When Should You Adjust Standard Errors for Clustering?

Alberto Abadie		Susan Athey
MIT		Stanford
Guido W. Imbens		Jeffrey M. Wooldridge
Stanford		MSU

March 14, 2024

Abstract

Clustered standard errors, with clusters defined by factors such as geography, are widespread in empirical research in economics and many other disciplines. Formally, clustered standard errors adjust for the correlations induced by sampling the outcome variable from a data-generating process with unobserved cluster-level components. However, the standard econometric framework for clustering leaves important questions unanswered: (i) Why do we adjust standard errors for clustering in some ways but not others, e.g., by state but not by gender, and in observational studies, but not in completely randomized experiments? (ii) Why is conventional clustering an “all-or-nothing” adjustment, while within-cluster correlations can be strong or extremely weak? (iii) In what settings does the choice of whether and how to cluster make a difference? We address these and other questions using a novel framework for clustered inference on average treatment effects. In addition to the common sampling component, the new framework incorporates a design component that accounts for the variability induced on the estimator by the treatment assignment mechanism. We show that, when the number of clusters in the sample is a non-negligible fraction of the number of clusters in the population, conventional cluster standard errors can be severely inflated, and propose new variance estimators that correct for this bias.

1 Introduction

Imagine you estimated the effect of attending college on labor earnings using linear regression on a cross-section of U.S. workers. How should you calculate the standard error? Empirical studies in economics often report heteroskedasticity-robust standard errors (henceforth “robust”) associated with the work by eicker1963, huber1967behavior, and white1980heteroskedasticity. A common alternative is to report cluster-robust standard errors (henceforth “cluster”) associated with the work by liang1986longitudinal and arellano1987practitioners, with clustering often applied within geographic units such as states or counties. moulton1986random, moulton1987diagnostics and Bertrand2004did have shown that clustering adjustments can make a substantial difference, and since the 1980s cluster standard errors have become commonplace in empirical economics.

Later in this section, we estimate a log-linear regression of earnings on an indicator for some college using data from the 2000 U.S. Census. We find that standard errors clustered at the state level are more than 20 times larger than robust standard errors. Which ones should a researcher report? The conventional framework for clustering [see cameron2015practitioner, mackinnon2021cluster, for recent reviews] suggests that if the clustering adjustment matters, in the sense that the cluster standard errors are substantially larger than the robust standard errors, one should use the cluster standard errors. In this article, we develop a new framework for cluster adjustments to standard errors that nests the conventional framework as a limiting case. The new framework suggests novel standard error formulas that can substantially improve over robust and cluster standard errors in settings like the earnings regression described above.

Our proposed clustering framework differs from the standard one in that it includes a design component that accounts for between-clusters variation in treatment assignments. We argue that the new design component is important because between-cluster variation in treatment assignments often motivates the use of clustered standard errors in empirical studies [see, e.g., gentzkow2008preschool, cohen2010free]. In addition, our framework shifts the focus of interest from features of infinite super-populations/data-generating processes to average treatment effects defined for the finite (but potentially large) population at hand. As a result of this shift, it is the sampling process and the treatment assignment mechanism that solely determine the correct level of clustering; the presence of cluster-level unobserved components of the outcome variable becomes irrelevant for the choice of clustering level. Moreover, by focusing on finite populations (which could be entirely or substantially sampled in the data) we obtain standard errors smaller than those aiming to measure uncertainty with respect to features of infinite super-populations. We derive the large sample variances for the least squares and fixed effect estimators under our proposed framework and show that they differ in general from both the robust and the cluster variances. We also propose two estimators for the large sample variances, one analytic and one based on a re-sampling (bootstrap) approach. For the U.S. earnings application, our proposals produce standard errors that are substantially larger than the robust standard errors, but also substantially smaller than the conventional version of cluster standard errors.

We use our framework to highlight three common misconceptions surrounding clustering adjustments. The first misconception is that the need for clustering hinges on the presence of a non-zero correlation between residuals for units belonging to the same cluster. We show that the presence of such correlation does not imply the need to use cluster adjustments, and that the absence of such correlation does not imply that clustering is not required. The second misconception is that there is no harm in using clustering adjustments when they are not required, with the implication that if clustering the standard errors makes a difference, one should do so. To see that both of these claims are in fact incorrect, consider the following simple example. Suppose that, based on a random sample from the population of interest, we use the sample average of a variable to estimate its population mean. Suppose also that the population can be partitioned into clusters such as geographical units. If outcomes are positively correlated within clusters, the cluster variance will be larger than the robust variance. However, standard sampling theory directly implies that if the units are sampled randomly from the population there is no need to cluster. The harm in clustering in this case is that confidence intervals will be unnecessarily conservative, possibly by a wide margin. A third misconception is that researchers have only two choices: either fully adjust for clustering and use the cluster standard errors, or not adjust the standard errors at all and use the robust standard errors. We show that a combination of the robust and the cluster variance estimators can substantially improve accuracy over its two components.

The new clustering framework in this article has the advantage of providing actionable guidance on a question of substantial consequence for empirical practice in econometrics: When should standard errors be clustered, and at what level? In the conventional model-based econometric framework, the researcher takes a stand on the error component structure of a model for the outcome variable. For example, suppose that, following moulton1986random, moulton1987diagnostics, the researcher posits a random effects model, with random effects at the state level. In this setting, a repeated sampling thought experiment entails that, for each sample, different values of the state random effects are drawn from their distributions. This model-based approach implies that if we are estimating a population mean using a sample average one needs to cluster the standard errors at the state level even if the sample is a random sample of individuals and not a clustered sample. A drawback of the model-based econometric framework for clustering is that empirical researchers need to take a stand on the structure of the error components of their models.

A second, closely related, framework for clustering that is often invoked in the econometrics literature is motivated by a sampling mechanism that in a first stage selects clusters at random from an infinite population, followed by a second stage of random sampling of units from the sampled clusters (or keeping all units in a cluster). Although this framework is appropriate for some applications in the analyses of surveys, where it originated [kish1995survey, thompson2012sampling], we argue that it is not appropriate for many of the data sets economists and other social scientists analyze. In many applications in economics, researchers do observe units from all the clusters they are interested in, e.g., all the states in the U.S., and a framework based on randomly sampling a small fraction of a large population of clusters does not apply.

Neither of the two conventional frameworks for clustered inference described above fully incorporates the design aspect of clustering. And it is the lack of a design component that makes them inappropriate for inference on treatment effects. To gain insight on the importance of the assignment mechanism for the standard errors of treatment effects estimators, consider a setting with individuals sampled at random from a population, but where treatment is assigned at the cluster level, with the same treatment value for all the individuals in the same cluster. Assume that the quantity of interest is the population average treatment effect. Clustered assignment to treatment is equivalent to clustered sampling of potential outcomes. Because the parameter of interest depends on averages of potential outcomes, which are sampled in a clustered manner, clustering of the standard errors is required in this setting, even when the individual observations are sampled at random. Our framework for clustered inference in this setting is close in spirit to the sampling framework described in the previous paragraph, but it incorporates a design component.

By shifting the attention from parameters of a data generating process for the outcomes to the average treatment effect for the population at hand, a researcher applying the proposals in this article does not need to take a stand on the error component structure of a model for the outcome variable to calculate standard errors. Instead, all the relevant variability of the estimator with respect to the average treatment effect is generated by the sampling mechanism, which extracts the sample from the population, and the assignment mechanism, which determines which units are exposed to the treatment. We see this as an intrinsic advantage of the framework proposed in this article in settings where it is difficult to justify a particular error component structure.

In this article we make three contributions. The first one is a novel framework for clustering, building on the one developed by abadie2020sampling for the analysis of regression estimators from a design perspective. We allow for clustering both in the sampling process and in the assignment process. As a result, the framework nests both the traditional case of clustered sampling and the case of clustered treatment assignment in experiments as special cases. It also allows for intermediate cases. In particular, treatment assignment may depend on cluster but not perfectly so, and there remains variation in treatments within-clusters. This framework clarifies the separate roles of clustering in the sampling process and clustering in the assignment process. It also clarifies what we can learn from the data about the need to adjust standard errors for clustering. In our framework, the data are not informative about the need to adjust for clustering in the sampling process, but they are informative about the need to adjust for clustering in the assignment process.

In our second contribution, we derive central limit theorems and large sample variances for the least squares and the fixed effect estimators of average treatment effects that take into account variation both from sampling and assignment. Comparing these variances to limit versions of the robust and cluster variances shows that the robust standard errors are generally too small, and the cluster standard errors are unnecessarily conservative. These comparisons also highlight how heterogeneity in treatment effects affects inference in the estimation of average treatment effects. Often researchers specify models that implicitly assume constant treatment effects without appreciating the implications for inference. We show, however, that heterogeneity in treatment effects introduces additional variance components that affect the need for clustering adjustments.

In our third contribution, we propose new variance formulas and bootstrap procedures for treatment effects estimators in the presence of clustering. We use the term Causal Cluster Variance (CCV) for the analytic variance formulas. For the case of a least squares estimator of average treatment effects, the intuition for the CCV variance formula is as follows. The error of the least squares estimator is approximately equal to a sum, over all units, of an expression involving products of regression errors and regressors values. The robust variance is approximately equal to a sum, over all units, of the squares of these products. In contrast, the conventional cluster variance estimator is approximately equal to a sum, over all clusters, of squares of within-cluster sums of the same products. Although the sum over all clusters of the expectation of the within-cluster sums of these products is zero, the expectation for each cluster separately is not. For each cluster in the sample, it is possible to estimate the expectation of the sum of the products between regression errors and regressors values. The CCV formula uses these estimates to correct the bias of the conventional cluster variance. The CCV correction does not help much if only a small fraction of clusters are sampled. However, when a large fraction of the clusters are represented in the sample, the CCV correction can lead to substantial improvements. This adjustment relies on estimates of cluster-level treatment effects, and thus requires within-cluster variation in treatment assignment. In addition, we propose a bootstrap version of the variance estimator., which we compare to two benchmarks. In contrast to conventional bootstrap procedures, which are based on resampling individual units or entire clusters of units, our proposed Two-Stage-Cluster-Bootstrap (TSCB) conducts resampling in two stages. In the first stage, the fraction treated for each cluster is drawn from the empirical distribution of cluster-specific treatment fractions. In the second stage, the researcher samples the treated and control units from each cluster, with their number of units determined in the first stage. The CCV and TSCB variance estimators are designed for applications with large number of observations and substantial variation in treatment assignment within clusters.

To illustrate the empirical relevance of our results, we analyze a sample from the 2000 U.S. Decennial Census, which includes 2,632,838 individuals. We define 52 clusters according to residency in the 50 states, Puerto Rico, and the District of Columbia. We consider two log-linear regressions of individual earnings on a treatment variable that encodes information on college attendance. In the first specification, the treatment variable is measured as an average, at the state level. In a second specification, we measure college attendance at the individual level.

Table 1: College effects in the Census sample

Dependent variable: Log labor earnings
Panel A
Treatment:	State indicator for share of some
	college greater than 0.55
		OLS
coefficient		0.1022
standard error:
robust		(0.0012)
cluster		(0.0312)
Panel B
Treatment:	Individual indicator for some college
		OLS		FE
coefficient		0.4656		0.4570
standard error:
robust		(0.0012)		(0.0012)
cluster		(0.0269)		(0.0276)
causal cluster variance (CCV)		(0.0035)		(0.0014)
two-stage cluster bootstrap (TSCB)		(0.0036)		(0.0014)

In Panel A of Table 1, we report results for a regression where the only explanatory variable is a binary treatment that takes value one if the fraction of individuals with at least some college residing in the state is 0.55 or higher, and value zero otherwise (we chose the 0.55 value to ensure sufficient variation in the treatment over the 52 clusters). Notice that the treatment is constant within states. We report the ordinary least squares (OLS) estimate, as well as robust and cluster standard errors. Since the late 1980s, it has been common practice to report cluster standard errors in settings where the regressors are constant within a cluster. Clustering at the state level makes a substantial difference relative to using robust standard errors, with the cluster standard errors approximately twenty-six times larger than the robust standard errors.

In Panel B of Table 1, the sole regressor is an individual-level indicator for at least some college. In addition to OLS, we report the fixed effects (FE) estimate (with fixed effects for the 50 states, plus Washington DC and Puerto Rico) and robust, cluster, CCV, and TSCB standard errors in parentheses. Like for the regression of the first panel, clustering at the state level makes a substantial difference in the standard errors, with the cluster standard errors approximately twenty-three times larger than the robust standard errors, both for the OLS and the FE regressions. In Panel B, our proposed CCV and TSCB standard errors for the OLS estimate are 0.0035 and 0.0036 respectively, in between the robust standard errors (0.0012) and the cluster standard errors (0.0269), and substantially different from both. The same holds for the FE estimator. The cluster standard error is 0.0276, quite different from the robust standard errors, 0.0012. The CCV and TSCB standard errors are 0.0014, in between robust and cluster but much closer to robust.

2 A Framework for Clustering

In this section, we describe in detail the framework for our analysis. There are multiple components to our set-up that are not explicitly modeled in the usual analysis of the variance of econometric estimators. In general, quantifying the uncertainty of parameter estimates requires describing the population and articulating the assumptions that describe how the sample was generated from that population (that is, building a model for the data generating process). In our framework, there are three distinct sources of sampling variation that lead to variation in the estimates. First, there is variation across samples in which units are observed in each cluster. Second, there is potentially variation in which clusters are observed (which leads to different units being observed). Third, there is variation in the treatment assignment across units. Whereas the standard framework for clustering focuses solely on the first two (sampling) sources of uncertainty, our proposed framework allows for all three. How much these three components matter for the variance of the least squares and fixed effects estimators of the average treatment effect depends on (i) the sampling process, (ii) the assignment process, and (iii) the heterogeneity in the treatment effects across clusters. To facilitate the calculation of asymptotic approximations in a range of relevant settings for empirical practice, it is convenient to formally consider a sequence of populations where we can separately control the fraction of units in the population that are sampled and the fraction of clusters in the population that is sampled, as well as the assignment mechanism.

2.1 A Sequence of Populations

We have a sequence of populations indexed by $k$. The $k$-th population has $n_{k}$ units, indexed by $i=1,\ldots,n_{k}$. The population is partitioned into $m_{k}$ clusters. Let $m_{k,i}\in\{1,\ldots,m_{k}\}$ denote the cluster that unit $i$ of population $k$ belongs to. The number of units in cluster $m$ of population $k$ is $n_{k,m}\geq 1$. For each unit, $i$, there are two potential outcomes, $y_{k,i}(1)$ and $y_{k,i}(0)$, corresponding to treatment and no treatment. Thus the population is characterized by the set of triples $(m_{k,i},y_{k,i}(0),y_{k,i}(1))$, for units $1,\ldots,n_{k}$ and clusters $1,\ldots,m_{k}$. The object of interest is the population average treatment effect

$$\tau_{k}=\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}\bigl{(}y_{k,i}(1)-y_{k,i}(0)\bigr{)}.$$

The population average treatment effect by cluster is

$$\tau_{k,m}=\frac{1}{n_{k,m}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}(y_{k,i}(1)-y_{k,i}(0)).$$

Therefore,

$$\tau_{k}=\sum_{m=1}^{m_{k}}\frac{n_{k,m}}{n_{k}}\tau_{k,m}.$$

We assume that potential outcomes, $y_{k,i}(1)$ and $y_{k,i}(0)$, are bounded in absolute value, uniformly for all $(k,i)$.

For each unit in the population, we define the stochastic treatment indicator, $W_{k,i}\in\{0,1\}$. The realized outcome for unit $i$ in population $k$ is $Y_{k,i}=y_{k,i}(W_{k,i})$. For a random sample of the population, we observe the triple $(Y_{k,i},W_{k,i},m_{k,i})$. Inclusion in the sample is represented by the random variable $R_{k,i}$, which takes value one if unit $i$ belongs to the sample, and value zero if not. We next describe the two components of the stochastic nature of the sample: the sampling process that determines the values of $R_{k,i}$, and the assignment process that determines the values of $W_{k,i}$.

2.2 The Sampling Process

The sampling process that determines the values of $R_{k,i}$ is independent of the potential outcomes and the assignments. It consists of two stages. First, clusters are sampled with cluster sampling probability $q_{k}\in(0,1]$. Second, units are sampled from the subpopulation consisting of all the sampled clusters, with unit sampling probability equal to $p_{k}\in(0,1]$. Both $q_{k}$ and $p_{k}$ may be equal to one, or close to zero. If $q_{k}=1$, we sample all clusters. If $p_{k}=1$, we sample all units from the sampled clusters. If $q_{k}=p_{k}=1$, all units in the population are sampled. The standard framework for analyzing clustering focuses on the special case where $q_{k}\rightarrow 0$, so only a small fraction of the clusters in the population are sampled. The case $q_{k}=1$ and $p_{k}\rightarrow 0$ corresponds to taking a relatively small random sample of units from the population. While this is an important special case, there are also many applications where the sampled clusters comprise a large fraction of the overall set of clusters. We refer to the case of $q_{k}=1$ as random sampling and to the case of $q_{k}<1$ as clustered sampling.

2.3 The Assignment Process

The assignment process that determines the values of $W_{k,i}$ also consists of two stages. In the first stage of the assignment process, for cluster $m$ in population $k$, an assignment probability $A_{k,m}\in[0,1]$ is drawn randomly from a distribution with mean $\mu_{k}$, bounded away from zero and one uniformly in $k$, and variance $\sigma^{2}_{k}$, independently for each cluster. The variance $\sigma^{2}_{k}$ is key. If $\sigma^{2}_{k}$ is zero, then $A_{k,m}$ is the same for all $m$, and $W_{k,i}$ is randomly assigned across clusters. We refer to this case as random assignment. For positive values of $\sigma_{k}^{2}$ assignment probabilities depend on cluster. Because $A_{k,m}^{2}\leq A_{k,m}$, it follows that $\sigma_{k}^{2}$ is bounded above by $\mu_{k}(1-\mu_{k})$ and that the bound is attained when $A_{k,m}$ can only take values zero or one, so all units within a cluster have the same values for the treatment. We use the term clustered assignment to refer to the case $\sigma_{k}^{2}=\mu_{k}(1-\mu_{k})$, when there is no within-cluster variation in $W_{k,i}$. We use the term partially clustered assignment to refer to the case $0<\sigma_{k}^{2}<\mu_{k}(1-\mu_{k})$, where assignment depends on cluster but not all units in the same cluster necessarily have the same value of $W_{k,i}$. In the second stage of the assignment process, each unit in cluster $m$ is assigned to the treatment independently, with cluster-specific probability $A_{k,m}$.

3 The Least Squares Estimator and its Variance

Let

$$N_{k,1}=\sum_{i=1}^{n_{k}}R_{k,i}W_{k,i}\quad\mbox{ and }\quad N_{k,0}=\sum_{i=1}^{n_{k}}R_{k,i}(1-W_{k,i})$$

be the number of treated and untreated units in the sample, respectively; these are random variables. The total sample size is $N_{k}=N_{k,1}+N_{k,0}$.

We first analyze the OLS estimator of a regression of the outcome $Y_{k,i}$ on an intercept and the treatment indicator $W_{k,i}$. The OLS estimator (modified so it is well-defined even when $N_{k,1}=0$ or $N_{k,0}=0$) is equal to the difference in means:

$$\widehat{\tau}_{k}=\frac{1}{N_{k,1}\vee 1}\sum_{i=1}^{n_{k}}R_{k,i}W_{k,i}Y_{k,i}-\frac{1}{N_{k,0}\vee 1}\sum_{i=1}^{n_{k}}R_{k,i}(1-W_{k,i})Y_{k,i},$$ (1)

where $N_{k,1}\vee 1$ and $N_{k,0}\vee 1$ are the maxima of $N_{k,1}$ and 1 and of $N_{k,0}$ and 1, respectively.

We make the following assumptions about the sampling process and the cluster sizes: (i) $m_{k}q_{k}\rightarrow\infty$, (ii) $\liminf_{k\rightarrow\infty}p_{k}\min_{m}n_{k,m}>0$, and (iii) $\limsup_{k\rightarrow\infty}\max_{m}n_{k,m}/\min_{m}n_{k,m}<\infty$. The first assumption implies that the expected number of sampled clusters goes to infinity as $k$ increases. The second assumption implies that the average number of observations sampled per cluster, conditional on the cluster being sampled, does not go to zero. The third assumption restricts the imbalance between the number of units across clusters. Notice that assumptions (i) and (ii) imply $n_{k}p_{k}q_{k}\rightarrow\infty$, so the sample size becomes larger in expectation as $k$ increases.

3.1 Large $k$ Distribution of the Least Squares Estimator

Our first main result derives the large $k$ distribution of $\widehat{\tau}_{k}$. Let $\alpha_{k}=(1/n_{k})\sum_{i=1}^{n_{k}}y_{k,i}(0)$, $u_{k,i}(1)=y_{k,i}(1)-(\alpha_{k}+\tau_{k})$, and $u_{k,i}(0)=y_{k,i}(0)-\alpha_{k}$. Under additional regularity conditions in the Appendix,

$$\sqrt{N_{k}}(\widehat{\tau}_{k}-\tau_{k})/v_{k}^{1/2}\stackrel{{\scriptstyle d}}{{\longrightarrow}}N(0,1),$$

where

$$v_{k} =\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}\bigg{(}\frac{u^{2}_{k,i}(1)}{\mu_{k}}+\frac{u^{2}_{k,i}(0)}{1-\mu_{k}}\bigg{)}$$

$$-p_{k}\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}\big{(}u_{k,i}(1)-u_{k,i}(0)\big{)}^{2}-p_{k}\sigma_{k}^{2}\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}\bigg{(}\frac{u_{k,i}(1)}{\mu_{k}}+\frac{u_{k,i}(0)}{1-\mu_{k}}\bigg{)}^{2}$$

$$+p_{k}(1-q_{k})\frac{1}{n_{k}}\sum_{m=1}^{m_{k}}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}\big{(}u_{k,i}(1)-u_{k,i}(0)\big{)}\Bigg{)}^{2}$$

$$+p_{k}\sigma^{2}_{k}\frac{1}{n_{k}}\sum_{m=1}^{m_{k}}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}\bigg{(}\frac{u_{k,i}(1)}{\mu_{k}}+\frac{u_{k,i}(0)}{1-\mu_{k}}\bigg{)}\Bigg{)}^{2}.$$ (2)

The expression for the variance $v_{k}$ has multiple terms that make its interpretation challenging. We first interpret $v_{k}$ in some special cases to highlight the implications of clustered sampling and clustered assignment. In Section 3.3, we compare $v_{k}$ to the large-$k$ form of the robust and cluster variance estimators.

For the case of random sampling ($q_{k}=1$) and random assignment ($\sigma_{k}^{2}=0$), the variance simplifies to

$$\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}\bigg{(}\frac{u^{2}_{k,i}(1)}{\mu_{k}}+\frac{u^{2}_{k,i}(0)}{1-\mu_{k}}\bigg{)}-p_{k}\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}\big{(}u_{k,i}(1)-u_{k,i}(0)\big{)}^{2}.$$

As we show in Section 3.2 below, the first term in this variance is estimated by the robust variance estimator. The second term is a finite sample correction that is familiar from the literature on randomized experiments [e.g., neyman1923, imbens2015causal, abadie2020sampling]. This finite sample correction vanishes if there is either no heterogeneity in the treatment effects (so $u_{k,i}(1)-u_{k,i}(0)=y_{k,i}(1)-y_{k,i}(0)-\tau_{k}=0$), or if the sample is a small fraction of the population ($p_{k}\approx 0$).

Adding clustered sampling, $q_{k}<1$, increases the variance by

$$p_{k}(1-q_{k})\frac{1}{n_{k}}\sum_{m=1}^{m_{k}}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}\big{(}u_{k,i}(1)-u_{k,i}(0)\big{)}\Bigg{)}^{2},$$

which is the same as

$$p_{k}(1-q_{k})\frac{1}{n_{k}}\sum_{m=1}^{m_{k}}n_{k,m}^{2}(\tau_{k,m}-\tau_{k})^{2}.$$

This term vanishes if there is no heterogeneity in the average treatment effect across clusters. Although the sample is informative about heterogeneity in cluster average treatment effects, it is not informative about the value of $q_{k}$. Information about the need to adjust for clustered sampling ($q_{k}<1$) must come from outside the sample.

Clustered assignment, $\sigma^{2}_{k}>0$, adds two terms to the variance,

$$-p_{k}\sigma_{k}^{2}\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}\bigg{(}\frac{u_{k,i}(1)}{\mu_{k}}+\frac{u_{k,i}(0)}{1-\mu_{k}}\bigg{)}^{2}\!+p_{k}\sigma^{2}_{k}\frac{1}{n_{k}}\sum_{m=1}^{m_{k}}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}\bigg{(}\frac{u_{k,i}(1)}{\mu_{k}}+\frac{u_{k,i}(0)}{1-\mu_{k}}\bigg{)}\Bigg{)}^{2}.$$

As we explain in more detail in section 3.3, the sign of this expression depends on the amount of variation in potential outcomes that can be explained by the clusters. Note that in contrast to the lack of sample information about the need to adjust for clustered sampling, the sample is potentially informative about the need to account for clustered assignment.

The five terms making up the asymptotic variance $v_{k}$ can be of different order. The first term is an average of bounded terms, and so under our assumptions will be of order $\mathcal{O}(1)$. The second and third terms will be at most of the same order as the first one. If $p_{k}\approx 0$ so we can think of the sample as small relative to the population of sampled clusters, the first term dominates the second and third terms. If cluster sizes are bounded as $k$ increases, the fourth and fifth terms in are also order $\mathcal{O}(1)$. If, on the other hand, cluster sizes increase with $k$, these terms can be of higher order and dominate the variance. Whether they do so or not depends on the (i) magnitude of $p_{k}$, (ii) presence of clustering in sampling, (iii) presence of clustering in assignment, and (iv) heterogeneity in potential outcomes.

3.2 The Robust and Cluster Robust Variance Estimators

Let $\widehat{U}_{k,i}=Y_{k,i}-\widehat{\alpha}_{k}-\widehat{\tau}_{k}W_{k,i}$ be the residuals from the regression of $Y_{k,i}$ or a constant and $W_{k,i}$. Here, $\widehat{\alpha}_{k}$ is the intercept of the regression and $\widehat{\tau}_{k}$ is the coefficient on $W_{k,i}$ (equal to the expression in (1) with probability approaching one).

There are two common estimators of the variance of $\sqrt{N_{k}}(\widehat{\tau}_{k}-\tau_{k})$. First, the conventional robust variance estimator (eicker1963, huber1967behavior, white1980heteroskedasticity):

$$\widehat{V}^{\rm robust}_{k}=\frac{1}{\hbox{\vbox{\hrule height=0.5pt\kern 1.29167pt\hbox{\kern-0.20004pt$W$\kern-0.20004pt}}}_{k}^{2}(1-\hbox{\vbox{\hrule height=0.5pt\kern 1.29167pt\hbox{\kern-0.20004pt$W$\kern-0.20004pt}}}_{k})^{2}}\left\{\frac{1}{N_{k}}\sum_{i=1}^{n_{k}}R_{k,i}\widehat{U}_{k,i}^{2}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.29167pt\hbox{\kern-0.20004pt$W$\kern-0.20004pt}}}_{k})^{2}\right\},$$ (3)

where

$$\hbox{\vbox{\hrule height=0.5pt\kern 1.29167pt\hbox{\kern-0.20004pt$W$\kern-0.20004pt}}}_{k}=\frac{1}{N_{k}\vee 1}\sum_{i=1}^{n_{k}}R_{k,i}W_{k,i}.$$

Let

$$v^{\rm robust}_{k}=\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}\bigg{(}\frac{u^{2}_{k,i}(1)}{\mu_{k}}+\frac{u^{2}_{k,i}(0)}{1-\mu_{k}}\bigg{)}.$$

Under regularity conditions (see appendix), $\widehat{V}_{k}^{\rm{robust}}$ and $v_{k}^{\rm robust}$ are close in the following sense,

$$\frac{\widehat{V}_{k}^{\rm{robust}}}{v_{k}}=\frac{v_{k}^{\rm robust}}{v_{k}}+\scaleto{\mathcal{O}}{5pt}_{p}(1),$$

motivating our focus on the comparison of $v_{k}^{\rm robust}$ and $v_{k}$. In general the difference $v^{\rm robust}_{k}-v_{k}$ can be positive or negative, so the robust variance estimator can be invalid in large samples.

The second common variance estimator is the cluster variance [liang1986longitudinal, arellano1987practitioners],

$$\widehat{V}^{\rm cluster}_{k}=\frac{1}{\hbox{\vbox{\hrule height=0.5pt\kern 1.29167pt\hbox{\kern-0.20004pt$W$\kern-0.20004pt}}}_{k}^{2}(1-\hbox{\vbox{\hrule height=0.5pt\kern 1.29167pt\hbox{\kern-0.20004pt$W$\kern-0.20004pt}}}_{k})^{2}}\left\{\frac{1}{N_{k}}\sum_{m=1}^{m_{k}}\left(\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}\widehat{U}_{k,i}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.29167pt\hbox{\kern-0.20004pt$W$\kern-0.20004pt}}}_{k})\right)^{2}\right\}.$$ (4)

Define

$$v^{\rm cluster}_{k} =\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}\bigg{(}\frac{u^{2}_{k,i}(1)}{\mu_{k}}+\frac{u^{2}_{k,i}(0)}{1-\mu_{k}}\bigg{)}$$

$$-p_{k}\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}\big{(}u_{k,i}(1)-u_{k,i}(0)\big{)}^{2}-p_{k}\sigma_{k}^{2}\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}\bigg{(}\frac{u_{k,i}(1)}{\mu_{k}}+\frac{u_{k,i}(0)}{1-\mu_{k}}\bigg{)}^{2}$$

$$+p_{k}\frac{1}{n_{k}}\sum_{m=1}^{m_{k}}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}\big{(}u_{k,i}(1)-u_{k,i}(0)\big{)}\Bigg{)}^{2}$$

$$+p_{k}\sigma^{2}_{k}\frac{1}{n_{k}}\sum_{m=1}^{m_{k}}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}\bigg{(}\frac{u_{k,i}(1)}{\mu_{k}}+\frac{u_{k,i}(0)}{1-\mu_{k}}\bigg{)}\Bigg{)}^{2}.$$

Then, $\widehat{V}_{k}^{\rm{cluster}}$ is close to $v_{k}^{\rm cluster}$ in the sense that

$$\frac{\widehat{V}_{k}^{\rm{cluster}}}{v_{k}}=\frac{v_{k}^{\rm cluster}}{v_{k}}+\scaleto{\mathcal{O}}{5pt}_{p}(1).$$

The difference $v^{\rm cluster}_{k}-v_{k}$ is always nonnegative. Therefore, for large $k$, the cluster variance estimator can be conservative but cannot underestimate the variance of $\widehat{\tau}_{k}$.

3.3 Discussion

From the formulas for $v_{k}$, $v_{k}^{\rm{robust}}$, and $v_{k}^{\rm{cluster}}$ it follows that if $p_{k}$ is small enough, then $v_{k}^{\rm{robust}}$ and $v_{k}^{\rm{cluster}}$ are approximately equal to $v_{k}$. In this case, clustered sampling and clustered assignment do not matter much because the probability that two sample units belong to the same cluster is small.

The difference $v_{k}^{\rm{robust}}-v_{k}$ depends on two terms. The first term,

$$p_{k}\frac{1}{n_{k}}\Bigg{[}\sum_{i=1}^{n_{k}}\big{(}u_{k,i}(1)-u_{k,i}(0)\big{)}^{2}-(1-q_{k})\sum_{m=1}^{m_{k}}n_{k,m}^{2}(\tau_{k,m}-\tau_{k})^{2}\Bigg{]},$$ (5)

is equal to zero when treatment effects are constant (in which case, $u_{k,i}(1)-u_{k,i}(0)=0$ for $i=1,\ldots,n_{k}$ and $\tau_{k,m}-\tau_{k}=0$ for all $m=1,\ldots,m_{k}$). If all clusters are sampled, so $q_{k}=1$, and treatment effects are heterogeneous, (5) is positive. When only a fraction of the clusters are sampled, $q_{k}<1$, the sign of (5) depends on the extent to which heterogeneity in treatment effects can be explained by the clusters. If there is no variation in average treatment effects across clusters, the expression in (5) is non-negative. However, when clusters explain much of the variation in treatment effects, the expression in (5) can be negative and very large in magnitude because of the factor $n_{k,m}^{2}$. The second term of $v_{k}^{\rm{robust}}-v_{k}$ is equal to

$$p_{k}\sigma^{2}_{k}\sum_{m=1}^{m_{k}}\frac{n_{k,m}}{n_{k}}\Bigg{[}\frac{1}{n_{k,m}}\sum_{i=1}^{n_{k}}1\{m_{k,i} =m\}\bigg{(}\frac{u_{k,i}(1)}{\mu_{k}}+\frac{u_{k,i}(0)}{1-\mu_{k}}\bigg{)}^{2}$$

$$-n_{k,m}\Bigg{(}\frac{1}{n_{k,m}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}\bigg{(}\frac{u_{k,i}(1)}{\mu_{k}}+\frac{u_{k,i}(0)}{1-\mu_{k}}\bigg{)}\Bigg{)}^{2}\Bigg{]}.$$ (6)

This term is equal to zero if there is no clustered assignment, that is, $\sigma^{2}_{k}=0$. If $\sigma^{2}_{k}>0$, the sign of (6) depends on how much of the heterogeneity in potential outcomes is explained by the clusters. The expression in (6) is close to zero when there is little heterogeneity in potential outcomes, so $u_{k,i}(1)$ and $u_{k,i}(0)$ are close to zero. If there is heterogeneity in potential outcomes but average potential outcomes are nearly constant across clusters, (6) is positive. When the clusters explain enough heterogeneity in potential outcomes (6) can be negative and potentially very large in magnitude because of the factor $n_{k,m}$ multiplying the second term of the sum in (6). That is, the robust variance formula can severely underestimate the variance of $\widehat{\tau}_{k}$.

Clustered standard errors are conservative in general, that is, $v_{k}^{\rm{cluster}}\geq v_{k}$. In particular, the difference $v_{k}^{\rm{cluster}}-v_{k}$ is

$$v_{k}^{\rm{cluster}}-v_{k}=p_{k}q_{k}\frac{1}{n_{k}}\sum_{m=1}^{m_{k}}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}\big{(}u_{k,i}(1)-u_{k,i}(0)\big{)}\Bigg{)}^{2},$$

which can be rewritten as

$$v_{k}^{\rm{cluster}}-v_{k}=\left(\frac{p_{k}n_{k}}{m_{k}}\right)q_{k}\left\{\frac{1}{m_{k}}\sum_{m=1}^{m_{k}}\left(\frac{n_{k,m}m_{k}}{n_{k}}\right)^{2}(\tau_{k,m}-\tau_{k})^{2}\right\}.$$ (7)

When the expected fraction of clusters in the sample, $q_{k}$, is small, or when the average treatment effect is nearly constant between clusters, then $v_{k}^{\rm{cluster}}\approx v_{k}$. Aside from these special cases, the $p_{k}n_{k}/m_{k}$ factor in the formula above indicates that cluster standard errors can be extremely conservative in general.

4 Two New Variance Estimators

Estimation of the variance of $\widehat{\tau}_{k}$ is challenging because the different terms in $v_{k}$ can be of different orders of magnitude. In this section, we propose two estimators of the variance of $\widehat{\tau}_{k}$ that allow us to correct the bias of the cluster variance estimator, one analytic, and one resampling-based. As the expression for the bias of the cluster variance in (7) shows, the cluster variance is heavily biased if the fraction of the sampled clusters is large and there is substantial variation in the cluster-specific treatment effects. Although the proposed analytic variance estimator is defined irrespective of the value of $\sigma_{k}^{2}$, in order to for the correction to be effective we need to be able to estimate the cluster-specific treatment effects, and thus we need $\sigma_{k}^{2}$ to be less than its maximum value of $\mu_{k}(1-\mu_{k})$ to ensure that there is variation in the treatment assignment within clusters. One of the proposed variance estimators is based on a correction to $\widehat{V}^{\rm cluster}_{k}$, and the other is based on resampling methods. An alternative would be to directly estimate the bias term in (7) and subtract that from the cluster variance. A challenge with this approach is that the estimation error for the adjustment term is large (often leading to negative variances estimates) because the order of magnitude of the correction is itself large and this approach did not work well in our simulations. We do not report formal results for the variance estimators in the current paper. We demonstrate their performance in the simulations in Section 6. There may well be further refinements possible.

If $q_{k}$ is close to zero, the proposed variance estimators are close to $\widehat{V}^{\rm cluster}_{k}$, which has little bias in that case. If $\sigma_{k}^{2}=\mu_{k}(1-\mu_{k})$ (that is, when $W_{k,i}$ is constant within clusters), the proposed resampling variance estimator is not defined. To be effective both variance estimators rely on estimating the variation in treatment effects across clusters, and therefore require a substantial number of both treated and control observations per cluster. The proposed variance estimators lead to substantial improvements over $\widehat{V}^{\rm cluster}_{k}$ in cases where $\widehat{V}^{\rm cluster}_{k}$ has a large upward bias. The downside of the proposed variance estimators is that they can be conservative when there is no need to cluster because there is no heterogeneity in treatment effects, or when there are too few treated and control observations per cluster to estimate the heterogeneity in the treatment effects precisely.

We first consider in Section 4.1 the case with $q_{k}=1$ so we have random sampling. Next we consider in Section 4.2 the case with clustered sampling $q_{k}<1$. In Section 4.3 we propose a bootstrap procedure for estimating the variance. The proposed variance estimators perform very well in the simulation study of Section 6. The derivation of their formal properties is left for future work.

4.1 The Case with All Clusters Observed

First we focus on the case with $q_{k}=1$ (all clusters observed), but allowing for general $p_{k}$. Let $U_{k,i}=W_{k,i}u_{k,i}(1)+(1-W_{k,i})u_{k,i}(0)$. The first step is to approximate the normalized error of the least squares estimator $\widehat{\tau}_{k}$ by a normalized sample average over clusters,

$$\sqrt{N_{k}}(\widehat{\tau}_{k}-\tau_{k})/v_{k}^{1/2}=\frac{1}{\sqrt{n_{k}p_{k}v_{k}}\mu_{k}(1-\mu_{k})}\sum_{m=1}^{m_{k}}C_{k,m}+\scaleto{\mathcal{O}}{5pt}_{p}(1),$$ (8)

where the terms

$$C_{k,m}=\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}(W_{k,i}-\mu_{k})U_{k,i}$$

are independent across clusters. In the appendix, we show

$$\widehat{V}_{k}^{\rm{cluster}}/v_{k}=\frac{1}{n_{k}p_{k}v_{k}}\left(\frac{1}{\mu_{k}(1-\mu_{k})}\right)^{2}\sum_{m=1}^{m_{k}}C_{k,m}^{2}+\scaleto{\mathcal{O}}{5pt}_{p}(1).$$ (9)

The expectation of $C_{m,k}$ is

$$E[C_{k,m}]=n_{k,m}p_{k}\mu_{k}(1-\mu_{k})(\tau_{k,m}-\tau_{k}),$$

with sum over clusters

$$\sum_{m=1}^{m_{k}}E[C_{k,m}]=p_{k}\mu_{k}(1-\mu_{k})\sum_{m=1}^{m_{k}}n_{k,m}(\tau_{k,m}-\tau_{k})=0.$$ (10)

That is, although the sum of the expectations of $C_{k,m}$ over clusters is equal to zero, these expectations are not equal to zero in general for each cluster separately. Because $\mbox{var}(C_{k,m})\leq E[C_{k,m}^{2}]$, the first term on the right-hand side of (9) is conservative on expectation relative to the variance of $\sqrt{N_{k}}(\widehat{\tau}_{k}-\tau_{k})/v_{k}^{1/2}$, which explains the conservativeness of $\widehat{V}_{k}^{\rm{cluster}}$.

Because of (10), we can replace the terms $C_{k,m}$ in (8) by $C_{k,m}-E[C_{k,m}]=C_{k,m,1}+C_{k,m,2}$, where

$$C_{k,m,1} =\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}(R_{k,i}-p_{k})(\tau_{k,m}-\tau_{k})\mu_{k}(1-\mu_{k}),$$

and

$$C_{k,m,2} =\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}\Bigl{(}(W_{k,i}-\mu_{k})U_{k,i}-(\tau_{k,m}-\tau_{k})\mu_{k}(1-\mu_{k})\Bigr{)}.$$

Therefore,

$$\sqrt{N_{k}}(\widehat{\tau}_{k}-\tau_{k})/v_{k}^{1/2}=\frac{1}{\sqrt{n_{k}p_{k}v_{k}}\mu_{k}(1-\mu_{k})}\left(\sum_{m=1}^{m_{k}}C_{k,m,1}+\sum_{m=1}^{m_{k}}C_{k,m,2}\right)+\scaleto{\mathcal{O}}{5pt}_{p}(1).$$ (11)

It can be shown that $C_{k,m,1}$ and $C_{k,m,2}$ have means equal to zero and are uncorrelated. In addition, $C_{k,m,1}$ and $C_{k,m,2}$ are uncorrelated across clusters. The variance of $\sum_{m=1}^{m_{k}}C_{k,m,1}/(\sqrt{n_{k}p_{k}}\mu_{k}(1-\mu_{k}))$ is

$$(1-p_{k})\sum_{m=1}^{m_{k}}\frac{n_{k,m}}{n_{k}}(\tau_{k,m}-\tau_{k})^{2}.$$

Let $\widehat{\tau}_{k,m}$ be difference between the sample average of the outcome for treated and nontreated units in cluster $m$. A direct estimator the variance of $\sum_{m=1}^{m_{k}}C_{k,m,2}$ is

$$\sum_{m=1}^{m_{k}}\left(\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}\Bigl{(}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.29167pt\hbox{\kern-0.20004pt$W$\kern-0.20004pt}}}_{k})\widehat{U}_{k,i}-(\widehat{\tau}_{k,m}-\widehat{\tau}_{k})\hbox{\vbox{\hrule height=0.5pt\kern 1.29167pt\hbox{\kern-0.20004pt$W$\kern-0.20004pt}}}_{k}(1-\hbox{\vbox{\hrule height=0.5pt\kern 1.29167pt\hbox{\kern-0.20004pt$W$\kern-0.20004pt}}}_{k})\Bigr{)}\right)^{2},$$ (12)

In practice, the estimator in (12) is biased from the correlations between the estimation errors of its components. We apply sampling splitting to address this bias. We first split the sample randomly into two subsamples. Let ${Z}_{k,i}\in\{0,1\}$ be the indicator that unit $i$ belongs to the second subsample, and let $\overline{Z}_{k}$ be the mean of $Z_{k,i}$. Using the subsample with ${Z}_{k,i}=0$, we obtain estimates $\widehat{\tau}_{k,m}^{\,*}$, $\widehat{\alpha}_{k}^{\,*}$, and $\widehat{\tau}_{k}^{\,*}$ of $\tau_{k,m}$, $\alpha_{k}$, and $\tau_{k}$, respectively. Next, for observations with $Z_{k,i}=1$, we calculate the residuals $\widehat{U}_{k,i}^{\,*}=Y_{k,i}-\widehat{\alpha}_{k}^{\,*}-\widehat{\tau}_{k}^{\,*}W_{k,i}$. Finally, we estimate the normalized variance for the case with $q_{k}=1$ as

$$\widehat{V}^{\rm CCV}_{k}(1) =\frac{1}{N_{k}\hbox{\vbox{\hrule height=0.5pt\kern 1.29167pt\hbox{\kern-0.20004pt$W$\kern-0.20004pt}}}_{k}^{2}(1-\hbox{\vbox{\hrule height=0.5pt\kern 1.29167pt\hbox{\kern-0.20004pt$W$\kern-0.20004pt}}}_{k})^{2}}\sum_{m=1}^{m_{k}}\Bigg{[}\frac{1}{\hbox{\vbox{\hrule height=0.5pt\kern 1.29167pt\hbox{\kern-0.20004pt$Z$\kern-0.20004pt}}}_{k}^{2}}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}Z_{k,i}\Big{(}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.29167pt\hbox{\kern-0.20004pt$W$\kern-0.20004pt}}}_{k})\widehat{U}_{k,i}^{\,*}$$

$$-(\widehat{\tau}_{k,m}^{\,*}-\widehat{\tau}_{k}^{\,*})\hbox{\vbox{\hrule height=0.5pt\kern 1.29167pt\hbox{\kern-0.20004pt$W$\kern-0.20004pt}}}_{k}(1-\hbox{\vbox{\hrule height=0.5pt\kern 1.29167pt\hbox{\kern-0.20004pt$W$\kern-0.20004pt}}}_{k})\Big{)}\Bigg{)}^{2}$$

$$-\frac{1-\hbox{\vbox{\hrule height=0.5pt\kern 1.29167pt\hbox{\kern-0.20004pt$Z$\kern-0.20004pt}}}_{k}}{\hbox{\vbox{\hrule height=0.5pt\kern 1.29167pt\hbox{\kern-0.20004pt$Z$\kern-0.20004pt}}}_{k}^{2}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}Z_{k,i}\Big{(}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.29167pt\hbox{\kern-0.20004pt$W$\kern-0.20004pt}}}_{k})\widehat{U}_{k,i}^{\,*}-(\widehat{\tau}_{k,m}^{\,*}-\widehat{\tau}_{k}^{\,*})\hbox{\vbox{\hrule height=0.5pt\kern 1.29167pt\hbox{\kern-0.20004pt$W$\kern-0.20004pt}}}_{k}(1-\hbox{\vbox{\hrule height=0.5pt\kern 1.29167pt\hbox{\kern-0.20004pt$W$\kern-0.20004pt}}}_{k})\Big{)}^{2}\Bigg{]}$$

$$+(1-p_{k})\sum_{m=1}^{m_{k}}\frac{\hbox{\vbox{\hrule height=0.5pt\kern 1.29167pt\hbox{\kern-0.20004pt$N$\kern-0.20004pt}}}_{k,m}}{N_{k}}(\widehat{\tau}_{k,m}-\widehat{\tau}_{k})^{2},$$ (13)

where $\hbox{\vbox{\hrule height=0.5pt\kern 1.29167pt\hbox{\kern-0.20004pt$N$\kern-0.20004pt}}}_{k,m}$ is the size of the sample in cluster $m$. For clusters with no variation in the treatment variable, we replace $\widehat{\tau}_{k,m}$ in (13) with $\widehat{\tau}_{k}$. For clusters with no variation in the treatment variable for a particular subsample, we replace $\widehat{\tau}_{k,m}^{*}$ in (13) with $\widehat{\tau}_{k}^{*}$. We derive the form of the CCV estimator in the appendix. To improve the precision of $\widehat{V}^{\rm CCV}_{k}(1)$, we re-estimate it multiple times with new sample splits (new values for $Z_{k,i}$) and then average the corresponding variance estimators. In our simulations of section 6, we re-estimate the variance estimator four times, and use sample splits with in expectation an equal number of units in each subsample, so $E[\overline{Z}_{k}]=1/2$.

4.2 The Case When Not All Clusters Are Sampled

To motivate the modification of the variance estimator $\widehat{V}_{k}^{\rm CCV}(1)$ for the $q_{k}<1$ case, notice that

$$v_{k}(q_{k})-v^{\rm cluster}_{k}=q_{k}\times(v_{k}(1)-v^{\rm cluster}_{k}),$$

where $v_{k}(q_{k})$ denotes the value of the true variance $v_{k}$ evaluated at $q_{k}$. That is, the variance for the general $q_{k}$ case is a convex combination of the true variance at $q_{k}=1$ and the cluster variance,

$$v_{k}(q_{k})=q_{k}\times v_{k}(1)+(1-q_{k})\times v^{\rm cluster}_{k}.$$

Let $\widehat{q}_{k}$ be the ratio between the number of sampled clusters and the total number of clusters in the population. The proposed variance estimator, $\widehat{V}^{\rm CCV}_{k}$, is a convex combination of $\widehat{V}^{\rm CCV}_{k}(1)$ and $\widehat{V}^{\rm cluster}_{k}$ with weights $\widehat{q}_{k}$ and $1-\widehat{q}_{k}$,

$$\widehat{V}^{\rm CCV}_{k}=\widehat{q}_{k}\times\widehat{V}^{\rm CCV}_{k}(1)+(1-\widehat{q}_{k})\times\widehat{V}^{\rm cluster}_{k}.$$ (14)

Computation of $\widehat{q}_{k}$ requires knowledge of $m_{k}$, the total number of clusters in the population.

4.3 A Bootstrap Variance Estimator

In the previous sections, we have discussed an analytic variance estimator. Here we suggest a resampling-based variance estimator, initially for the case with $q_{k}=1$. Like the causal bootstrap in imbens2021causal, the proposed bootstrap procedure takes into account the causal nature of the estimand and creates bootstrap samples where units (in this case clusters) have different assignments and assignment probabilities than they have in the original sample. It differs from earlier bootstrap variance estimators for clustered settings [e.g., cameron2015practitioner, menzel2021bootstrap] in that it allows for the possibility that a large fraction of clusters are observed.

The specific resampling procedure, which we call the two-stage-cluster-bootstrap (TSCB), consists of two stages. For each of the clusters, let $\hbox{\vbox{\hrule height=0.5pt\kern 1.29167pt\hbox{\kern-0.20004pt$N$\kern-0.20004pt}}}_{k,m}$ be the cluster-level sample size and $\hbox{\vbox{\hrule height=0.5pt\kern 1.29167pt\hbox{\kern-0.20004pt$W$\kern-0.20004pt}}}_{k,m}=N_{k,m,1}/(\hbox{\vbox{\hrule height=0.5pt\kern 1.29167pt\hbox{\kern-0.20004pt$N$\kern-0.20004pt}}}_{k,m}\vee 1)$ the cluster-level fraction of treated units. In the first stage of the bootstrap procedure, for each cluster we draw ${\overline{W}}^{\,b}_{k,m}$ with replacement from the empirical distribution of the cluster-level fractions of treated units, that is with probability $1/m_{k}$ from the set $\{{\overline{W}}_{k,1},\ldots,{\overline{W}}_{k,m_{k}}\}$. In the second stage, we draw $\hbox{\vbox{\hrule height=0.5pt\kern 1.29167pt\hbox{\kern-0.20004pt$N$\kern-0.20004pt}}}_{k,m}\overline{W}^{\,b}_{k,m}$ units with replacement from the set of treated units in cluster $m$ and $\hbox{\vbox{\hrule height=0.5pt\kern 1.29167pt\hbox{\kern-0.20004pt$N$\kern-0.20004pt}}}_{k,m}(1-\overline{W}^{\,b}_{k,m})$ units with replacement from the set of untreated units in cluster $m$. In order for this to be well-defined we do need all the ${\overline{W}}_{k,1}$ to be strictly between zero and one. We do this for all clusters to create the bootstrap sample, and calculate the bootstrap standard errors as the standard deviation of the treatment effect estimates across bootstrap iterations.

Next, consider the case with $q_{k}<1$. In this case, we need to take into account the fact that we see a fraction of the clusters in the population. We follow the approach proposed in chao1985bootstrap. Suppose $q=1/2,$ so we observe half the clusters in the population. The bootstrap procedure first creates a pseudo population consisting of the original population of clusters, plus one additional replica of each cluster. Then, to get a bootstrap sample, we sample randomly, without replacement, from the clusters in this pseudo population. Given the clusters in the bootstrap sample, we proceed as before, and ultimately calculate the bootstrap variance as the variance of the estimator over the bootstrap samples. chao1985bootstrap provide details and extensions to the case for the case where $1/q_{k}$ is not an integer.

The algorithm for the TSCB is summarized here.

Algorithm 1 Two Stage Cluster Bootstrap

Input:

      Sample $(Y_{k,i},W_{k,i},m_{k,i})$

      Fraction sampled clusters $q_{k}$

      Number of bootstrap replications $B$

Stage 1:

1a: Create pseudo population by replicating each cluster $1/q_{k}$ times

1b: For each cluster in the pseudo population, calculate the assignment probability $\overline{W}_{k,m}$

1c: Create a bootstrap sample of clusters by randomly drawing clusters from the pseudo population from Stage 1a

1d: For each sampled cluster, draw an assignment probability $A_{k,m}$ from the empirical distribution of the $\overline{W}_{k,m}$ from Stage 1b

Stage 2:

2a: Randomly draw from the set of treated units in cluster $m$, $\lfloor N_{k,m}A_{k,m}\rfloor$ units

2b: Randomly draw from the set of control units in cluster $m$, $\lfloor N_{k,m}(1-A_{k,m})\rfloor$ units

Calculations:

For the units in the bootstrap sample constructed in Stage 2, collect the values for $(Y_{k,i},W_{k,i},m_{k,i})$ and calculate the least squares or fixed effect estimator

Calculate the standard deviation of the least squares or fixed effect estimator over the $B$ bootstrap samples

5 The Fixed Effect Estimator

In this section, we report results for the fixed effect estimator often used in empirical research in economics. arellano1987practitioners, Bertrand2004did, cameron2015practitioner and mackinnon2021cluster have pointed out that cluster adjustments may still be necessary in fixed effects regressions. However, a view of clustering based on models with cluster-specific variance components creates ambiguity in the role of clustered standard errors for estimators with cluster fixed effects, which are specifically aimed to absorb cluster-level variation.

We first characterize the fixed effect estimator and derive its large $k$ distribution. Then, we discuss the properties of the two conventional variance estimators, the robust and cluster robust variance estimators. As in the least squares case, we find that the robust standard errors may be too small and the cluster standard errors may be unnecessarily large, especially in cases when the number of observations per cluster is large. We propose CCV and TSCB variance estimators. The CCV estimator for fixed effects has a different form than the one for least squares in section A.4.

The fixed effect estimator is based on a regression of the outcome on the treatment indicator and indicators for each of the clusters in the sample. It can be written as the least squares estimate for a regression of the outcome on the treatment, with both variables measured in deviation from cluster means,

$$\widehat{\tau}_{k}^{\rm{\,fixed}} =\frac{\displaystyle\sum_{m=1}^{m_{k}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}Y_{k,i}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.29167pt\hbox{\kern-0.20004pt$W$\kern-0.20004pt}}}_{k,m})}{\displaystyle\sum_{m=1}^{m_{k}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.29167pt\hbox{\kern-0.20004pt$W$\kern-0.20004pt}}}_{k,m})}.$$ (15)

Like in section 3, we assume that that potential outcomes are bounded, $m_{k}q_{k}\rightarrow\infty$, and $\limsup_{k\rightarrow\infty}\max_{m}n_{k,m}/\min_{m}n_{k,m}<\infty$. In addition, we assume (i) $(m_{k}q_{k})/((p_{k}n_{k})/m_{k})\allowbreak\rightarrow 0$, and (ii) the supports of the cluster probabilities, $A_{k,m}$, are bounded away from zero and one (uniformly in $k$ and $m$). Assumption (i) restricts the focus of our analysis in this section to settings where the expected number of sampled clusters is small relative to the expected number of sampled observations per sampled cluster. Together with the previous assumptions, assumption (i) implies $(p_{k}n_{k})/m_{k}\rightarrow\infty$, $n_{k}p_{k}q_{k}\rightarrow\infty$, and $p_{k}\min_{m}n_{k,m}\rightarrow\infty$. This last result, along with assumption (ii), ensures that $\widehat{\tau}_{k}^{\rm{\,fixed}}$ in (15) is well-defined with probability approaching one.

Let $\alpha_{k,m}=(1/n_{k,m})\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}y_{k,i}(0)$. For an observation, $i$, with $m_{k,i}=m$, we define the within-cluster residuals $e_{k,i}(0)=y_{k,i}(0)-\alpha_{k,m}$ and $e_{k,i}(1)=y_{k,i}(1)-\tau_{k,m}-\alpha_{k,m}$. Let

$$\tilde{v}_{k}=f_{k}/(\mu_{k}(1-\mu_{k})-\sigma^{2}_{k})^{2}$$ (16)

where

$$f_{k} =E[A_{k,m}(1-A_{k,m})^{2}]\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}e_{k,i}^{2}(1)+E[A^{2}_{k,m}(1-A_{k,m})]\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}e_{k,i}^{2}(0)$$

$$-p_{k}E[A^{2}_{k,m}(1-A_{k,m})^{2}]\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}(e_{k,i}(1)-e_{k,i}(0))^{2}$$

$$+\Big{(}E[A_{k,m}(1-A_{k,m})]-(5+p_{k})E[A^{2}_{k,m}(1-A_{k,m})^{2}]$$

$$\qquad\qquad+2q_{k}(E[A_{k,m}(1-A_{k,m})])^{2}\Big{)}\sum_{m=1}^{m_{k}}\frac{n_{k,m}}{n_{k}}(\tau_{k,m}-\tau_{k})^{2}$$

$$+\Big{(}p_{k}E[A^{2}_{k,m}(1-A_{k,m})^{2}]-p_{k}q_{k}(E[A_{k,m}(1-A_{k,m})])^{2}\Big{)}\sum_{m=1}^{m_{k}}\frac{n^{2}_{k,m}}{n_{k}}(\tau_{k,m}-\tau_{k})^{2}.$$

Under additional regularity conditions, which are described in the Appendix, we obtain the large $k$ distribution of the fixed effects estimator,

$$\sqrt{N_{k}}(\widehat{\tau}_{k}^{\rm{\,fixed}}-\tau_{k})/\tilde{v}_{k}^{1/2}\stackrel{{\scriptstyle d}}{{\longrightarrow}}N(0,1).$$ (17)

Let $\widetilde{U}_{k,i}=\widetilde{Y}_{k,i}-\widehat{\tau}_{k}^{\,{\rm fixed}}\widetilde{W}_{k,i}$, where $\widetilde{Y}_{k,i}=Y_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.29167pt\hbox{\kern-0.20004pt$Y$\kern-0.20004pt}}}_{k,m_{k,i}}$, $\widetilde{W}_{k,i}=(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.29167pt\hbox{\kern-0.20004pt$W$\kern-0.20004pt}}}_{k,m_{k,i}})$. The robust estimator of the variance of $\sqrt{N_{k}}(\widehat{\tau}^{\,{\rm fixed}}_{k}-\tau_{k})$ is

$$\widetilde{V}_{k}^{\rm{robust}}=\left.\frac{1}{N_{k}}\sum_{i=1}^{n_{k}}R_{k,i}\widetilde{W}_{k,i}^{2}\widetilde{U}_{k,i}^{2}\right/\left(\frac{1}{N_{k}}\sum_{i=1}^{n_{k}}R_{k,i}\widetilde{W}_{k,i}^{2}\right)^{2}.$$ (18)

Now let,

$$\tilde{v}_{k}^{\rm robust}=f_{k}^{\rm robust}/(\mu_{k}(1-\mu_{k})-\sigma^{2}_{k})^{2}.$$

with

$$f_{k}^{\rm robust} =E[A_{k,m}(1-A_{k,m})^{2}]\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}e^{2}_{k,i}(1)+E[A^{2}_{k,m}(1-A_{k,m})]\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}e^{2}_{k,i}(0)$$

$$+E[A_{k,m}(1-A_{k,m})(1-3A_{k,m}(1-A_{k,m}))]\sum_{m=1}^{m_{k}}\frac{n_{k,m}}{n_{k}}(\tau_{k,m}-\tau_{k})^{2}.$$

Notice that all terms of $f_{k}^{\rm robust}$ are bounded. In the appendix, we show that

$$\widetilde{V}_{k}^{\rm{robust}}=\tilde{v}_{k}^{\rm robust}+\scaleto{\mathcal{O}}{5pt}_{p}(1).$$

The cluster variance estimator for fixed effects is

$$\widetilde{V}_{k}^{\rm{cluster}}=\left.\frac{1}{N_{k}}\sum_{m=1}^{m_{k}}\left(\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}\widetilde{W}_{k,i}\widetilde{U}_{k,i}\right)^{2}\right/\left(\frac{1}{N_{k}}\sum_{i=1}^{n_{k}}R_{k,i}\widetilde{W}_{k,i}^{2}\right)^{2}.$$ (19)

Let,

$$\tilde{v}_{k}^{\rm cluster}=f_{k}^{\rm cluster}/(\mu_{k}(1-\mu_{k})-\sigma^{2}_{k})^{2}.$$

with

$$f_{k}^{\rm cluster} =E[A_{k,m}(1-A_{k,m})^{2}]\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}e_{k,i}^{2}(1)+E[A^{2}_{k,m}(1-A_{k,m})]\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}e_{k,i}^{2}(0)$$

$$-p_{k}E[A^{2}_{k,m}(1-A_{k,m})^{2}]\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}(e_{k,i}(1)-e_{k,i}(0))^{2}$$

$$+(E[A_{k,m}(1-A_{k,m})]-(5+p_{k})E[A^{2}_{k,m}(1-A_{k,m})^{2}])\sum_{m=1}^{m_{k}}\frac{n_{k,m}}{n_{k}}(\tau_{k,m}-\tau_{k})^{2}$$

$$+p_{k}E[A^{2}_{k,m}(1-A_{k,m})^{2}]\sum_{m=1}^{m_{k}}\frac{n^{2}_{k,m}}{n_{k}}(\tau_{k,m}-\tau_{k})^{2}.$$

We obtain in the appendix,

$$\frac{\widetilde{V}_{k}^{\rm{cluster}}}{\tilde{v}_{k}}=\frac{\tilde{v}_{k}^{\rm cluster}}{\tilde{v}_{k}}+\scaleto{\mathcal{O}}{5pt}_{p}(1).$$

Similar to the least squares case, the robust variance can underestimate the true variance, and the cluster variance is generally too large. Our proposed variance estimator is a convex combination of $\widetilde{V}^{\rm cluster}_{k}$ and and $\widetilde{V}^{\rm robust}_{k}$, with the weights selected to correct the bias of the cluster variance estimator as $k$ increases (see appendix for details).

$$\widetilde{V}_{k}^{\rm CCV}=\widehat{\lambda}_{k}\widetilde{V}^{\rm cluster}_{k}+(1-\widehat{\lambda}_{k}){\widetilde{V}}^{\rm robust}_{k}.$$ (20)

where the estimated weight for the cluster variance is

$$\widehat{\lambda}_{k}=1-q_{k}\,\frac{\left(\displaystyle\frac{1}{M_{k}}\sum_{m=1}^{m_{k}}Q_{k,m}\hbox{\vbox{\hrule height=0.5pt\kern 1.29167pt\hbox{\kern-0.20004pt$W$\kern-0.20004pt}}}_{k,m}(1-\hbox{\vbox{\hrule height=0.5pt\kern 1.29167pt\hbox{\kern-0.20004pt$W$\kern-0.20004pt}}}_{k,m})\right)^{2}}{\displaystyle\frac{1}{M_{k}}\sum_{m=1}^{m_{k}}Q_{k,m}\hbox{\vbox{\hrule height=0.5pt\kern 1.29167pt\hbox{\kern-0.20004pt$W$\kern-0.20004pt}}}_{k,m}^{2}(1-\hbox{\vbox{\hrule height=0.5pt\kern 1.29167pt\hbox{\kern-0.20004pt$W$\kern-0.20004pt}}}_{k,m})^{2}},$$

where $Q_{k,m}$ is an indicator that takes value one if cluster $m$ of population $k$ is sampled, and $M_{k}=\sum_{m=1}^{m_{k}}Q_{k,m}$ is the total number of sampled clusters. The second factor in the second term approximately (that is, ignoring the variance of $\overline{W}_{k,m}$ conditional on $A_{k,m}]$) estimates the variance of $A_{k,m}(1-A_{k,m})$ divided by its second moment, so that

$$\tilde{\lambda}\approx 1-q_{k}\frac{V(A_{k,m}(1-A_{k,m}))}{E[(A_{k,m}(1-A_{k,m}))^{2}]}.$$

If there is no variation in $W_{k,i}$ within any of the clusters the fixed effect estimator is not defined, and neither is this variance estimator. In all other cases the variance estimator is well-defined.

We also consider a bootstrap standard error, based on the same resampling procedure described in Section 4.3.

6 Simulations

We next report simulation results that illustrate the performance of the proposed variance estimators relative to existing alternatives. To operate in an empirically relevant setting, we create an artificial population based on the Census data briefly described in the introduction, which contains information on log earnings, an indicator for college attendance, and an indicator for state of residence for 2,632,838 individuals.

For each individual in this population of 2,632,838 individuals, we define $m_{k,i}$ using state of residence (plus Washington, DC, and Puerto Rico), for a total of 52 clusters. We assign potential outcomes as $y_{k,i}(0)=Y_{k,i}-\widehat{\tau}_{k,m}W_{k,i}$ and $y_{k,i}(1)=Y_{k,i}+\widehat{\tau}_{k,m}(1-W_{k,i})$, so treatment effects are constant within clusters. We then repeatedly create samples from this population. Creating a sample requires fixing $p_{k}$, $q_{k}$, and fixing the distribution of $A_{k,m}$ and then drawing from the implied distribution for $R_{k,i}$ and $W_{k,i}$ to generate outcomes for all sampled units. In the baseline design, we set $p_{k}=q_{k}=1$, so we sample all $m_{k}=52$ clusters and all $n_{k}=2{,}632{,}838$ individuals in the population. For the assignment mechanism in the baseline design, we convert cluster means of the treatment variable into log-odds, $\widehat{\ell}_{k,m}=\ln(\hbox{\vbox{\hrule height=0.5pt\kern 1.29167pt\hbox{\kern-0.20004pt$W$\kern-0.20004pt}}}_{k,m}/(1-\hbox{\vbox{\hrule height=0.5pt\kern 1.29167pt\hbox{\kern-0.20004pt$W$\kern-0.20004pt}}}_{k,m}))$. Let $(\widehat{\mu}_{\ell},\widehat{\sigma}_{\ell})$ be the average and the sample standard deviation of $\widehat{\ell}_{k,m}$. We then draw $\ln(A_{k,m}/(1-A_{k,m}))$ for cluster $m$ from a normal distribution with expected value $\widehat{\mu}_{\ell}$ and standard deviation $\widehat{\sigma}_{\ell}$. Given the cluster assignment probability $A_{k,m}$, we assign the treatment in cluster $m$ by drawing from a binomial distribution with parameter $A_{k,m}$.

We calculate the standard deviation of the least squares and fixed effect estimators, normalized by the square root of the sample size, $N_{k}^{1/2}\mbox{s.d.}$, across 10,000 samples drawn according to the procedure outlined above. This is the benchmark against which we compare the various estimates of standard errors. For the least squares and the fixed effects estimators, respectively, we first calculate the (infeasible) asymptotic standard errors $v_{k}^{1/2}$ and $\widetilde{v}_{k}^{1/2}$ to benchmark the performance of the feasible variance estimators. Next, we calculate the averages across 10,000 simulations of the robust, cluster, CCV, and TCSB standard errors, where we use 100 bootstrap replications in each simulation. Table 6 reports the results. Table LABEL:table:coverage_rates reports coverage rates for 95 percent confidence intervals. In the design column of the two tables $\sigma_{\tau_{k}}$ is the standard deviation of the cluster average treatment effect.

Table 2: Average standard errors across simulations

					normalized standard error
		​$N_{k}^{1/2}\mbox{s.d.}$​	​$v_{k}^{1/2}$​	​$\widetilde{v}_{k}^{1/2}$​	​robust​	​cluster​	​CCV​	​TSCB​
Baseline design: $p_{k}=1$, $q_{k}=1$, $\sigma_{\tau_{k}}=.120$, $\sigma_{k}=.057$
	OLS​	5.91	5.90		1.90	44.86	6.32	5.80
	FE​	2.34		2.32	1.90	44.63	2.31	2.29
Second Design: $p_{k}=.1$, $q_{k}=1$, $\sigma_{\tau_{k}}=.120$, $\sigma_{k}=.057$
	OLS​	2.61	2.59		1.90	14.28	3.78	2.60
	FE​	1.95		1.95	1.90	14.21	1.95	1.94
Third Design: $p_{k}=.1$, $q_{k}=1$, $\sigma_{\tau_{k}}=.480$, $\sigma_{k}=.206$
	OLS​	14.50	14.17		1.98	56.46	13.70	14.33
	FE​	12.14		11.89	2.13	56.79	11.61	12.07
Fourth design: $p_{k}=.1$, $q_{k}=1$, $\sigma_{\tau_{k}}=0$, $\sigma_{k}=.206$
	OLS​	9.39	9.39		1.90	8.20	9.19	9.37
	FE​	2.04		2.04	2.04	1.97	2.04	2.09
Fifth design: $p_{k}=.1$, $q_{k}=1$, $\sigma_{\tau_{k}}=.480$, $\sigma_{k}=0$
	OLS​	1.95	1.97		1.97	56.42	4.53	2.04
	FE​	1.91		1.94	1.94	56.42	1.96	1.90

• ​​Notes: $N_{k}^{1/2}$s.d. is the standard deviation of the estimators over the simulations, multiplied by the square root of the sample size. $v_{k}^{1/2}$ is the square root of the asymptotic variance in equation (3.1). $\tilde{v}_{k}^{1/2}$ is the square root of the asymptotic variance of the fixed effect estimator in (16). The remaining four columns report average values of robust, cluster, CCV, and TSCB standard errors across simulations (multiplied by $N_{k}^{1/2}$). $p_{k}$ and $q_{k}$ are the unit and cluster sampling probabilities, respectively. $\sigma_{\tau_{k}}$ is the standard deviation of the cluster average treatment effect. $\sigma_{k}$ is the standard deviation across clusters of the treatment assignment probabilities.

For the baseline design, the normalized standard deviation of the least squares estimator is 5.91. This is well approximated by the asymptotic standard error, 5.90. The robust standard error is on average over the simulations 1.90, less than one-third of the normalized standard deviation of the estimator. The cluster standard error is far too large, on average 44.86, more than seven times the value of the normalized standard deviation. CCV improves considerably over robust and cluster. The average CCV standard error is 6.32, about 7 percent higher than the normalized standard deviation. The TSCB standard error is the most accurate, on average equal to 5.80. For the fixed effect estimator, the asymptotic standard error is again accurate. The robust standard error is about 16 percent too small, leading to a coverage rate for the nominal 95 percent confidence interval of 0.89 in Table LABEL:table:coverage_rates. The cluster standard error is too large by a factor of 20. CCV and TSCB standard errors closely approximate the normalized standard error.

When Should You Adjust Standard Errors for Clustering? Alberto Abadie, Susan Athey, Guido W. Imbens, and Jeffrey M. Wooldridge Current version: appendix

A.1 Setting and notation

We have a sequence of populations indexed by $k$. The $k$-th population has $n_{k}$ units, indexed by $i=1,\ldots,n_{k}$. The population is partitioned into $m_{k}$ strata or clusters. Let $m_{k,i}\in\{1,\ldots,m_{k}\}$ denote the stratum that unit $i$ of population $k$ belongs to. The number of units in cluster $m$ of population $k$ is $n_{k,m}\geq 1$. For each unit, $i$, there are two potential outcomes, $y_{k,i}(1)$ and $y_{k,i}(0)$, corresponding to treatment and no treatment. The parameter of interest is the population average treatment effect

$$\tau_{k}=\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}(y_{k,i}(1)-y_{k,i}(0)).$$

The population treatment effect by cluster is

$$\tau_{k,m}=\frac{1}{n_{k,m}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}(y_{k,i}(1)-y_{k,i}(0)).$$

Therefore,

$$\tau_{k}=\sum_{m=1}^{m_{k}}\frac{n_{k,m}}{n_{k}}\tau_{k,m}.$$

We will assume that potential outcomes, $y_{k,i}(1)$ and $y_{k,i}(0)$, are bounded in absolute value, uniformly for all $(k,i)$.

We next describe the two components of the stochastic nature of the sample. There is a stochastic binary treatment for each unit in each population, $W_{k,i}\in\{0,1\}$. The realized outcome for unit $i$ in population $k$ is $Y_{k,i}=y_{k,i}(W_{k,i})$. For a random sample of the population, we observe the triple $(Y_{k,i},W_{k,i},m_{k,i})$. Inclusion in the sample is represented by the random variable $R_{k,i}$, which takes value one if unit $i$ belongs to the sample, and value zero if not.

The sampling process that determines the values of $R_{k,i}$ is independent of the potential outcomes and the assignments. It consists of two stages. First, clusters are sampled with cluster sampling probability $q_{k}\in(0,1]$. Second, units are sampled from the subpopulation consisting of all the sampled clusters, with unit sampling probability equal to $p_{k}\in(0,1]$. Both $q_{k}$ and $p_{k}$ may be equal to one, or close to zero. If $q_{k}=1$, we sample all clusters. If $p_{k}=1$, we sample all units from the sampled clusters. If $q_{k}=p_{k}=1$, all units in the population are sampled.

The assignment process that determines the values of $W_{k,i}$ also consists of two stages. In the first stage, for cluster $m$ in population $k$, an assignment probability $A_{k,m}\in[0,1]$ is drawn randomly from a distribution with mean $\mu_{k}$, bounded away from zero and one uniformly in $k$, and variance $\sigma^{2}_{k}$, independently for each cluster. The variance $\sigma^{2}_{k}$ is key. If it is zero, we have random assignment across clusters. For positive values of $\sigma_{k}^{2}$ we have correlated assignment within the clusters. Because $A_{k,m}^{2}\leq A_{k,m}$, it follows that $\sigma_{k}^{2}$ is bounded above by $\mu_{k}(1-\mu_{k})$ and that the bound is attained when $A_{k,m}$ can only take values zero or one (so all units within a cluster have the same values for the treatment). In the second stage, each unit in cluster $m$ is assigned to the treatment independently, with cluster-specific probability $A_{k,m}$.

A.2 Base case: Difference in means

Let

$$N_{k,1}=\sum_{i=1}^{n_{k}}R_{k,i}W_{k,i}\quad\mbox{ and }\quad N_{k,0}=\sum_{i=1}^{n_{k}}R_{k,i}(1-W_{k,i})$$

be the number of treated and untreated units in the sample, respectively. The total sample size is $N_{k}=N_{k,1}+N_{k,0}$. We consider the simple difference of means between treated and non-treated, which is obtained as the coefficient on the treatment indicator in a regression of the outcome on a constant and the treatment,

$$\widehat{\tau}_{k}=\frac{1}{N_{k,1}\vee 1}\sum_{i=1}^{n_{k}}R_{k,i}W_{k,i}Y_{k,i}-\frac{1}{N_{k,0}\vee 1}\sum_{i=1}^{n_{k}}R_{k,i}(1-W_{k,i})Y_{k,i}.$$

We make the following assumptions about the sampling process and the cluster sizes: (i) $q_{k}m_{k}\rightarrow\infty$, (ii) $\liminf_{k\rightarrow\infty}p_{k}\min_{m}n_{k,m}>0$, and (iii) $\limsup_{k\rightarrow\infty}\max_{m}n_{k,m}/\min_{m}n_{k,m}<\infty$. The first assumption implies that the expected number of sampled clusters goes to infinity as $k$ increases. The second assumption implies that the average number of observations sampled per cluster, conditional on the cluster being sampled, does not go to zero. The third assumption restricts the imbalance between the number of units across clusters. Notice that assumptions (i) and (ii) imply $n_{k}p_{k}q_{k}\rightarrow\infty$, so the sample size becomes larger in expectation as $k$ increases.

A.2.1 Large $k$ distribution

Let $\alpha_{k}=(1/n_{k})\sum_{i=1}^{n_{k}}y_{k,i}(0)$ and $\tau_{k}=(1/n_{k})\sum_{i=1}^{n_{k}}(y_{k,i}(1)-y_{i,k}(0))$, $u_{k,i}(1)=y_{k,i}(1)-(\alpha_{k}+\tau_{k})$, and $u_{k,i}(0)=y_{k,i}(0)-\alpha_{k}$. Notice that,

$$\sum_{i=1}^{n_{k}}u_{k,i}(1)=\sum_{i=1}^{n_{k}}u_{k,i}(0)=0.$$

This implies

$$\sqrt{n_{k}p_{k}q_{k}}(\widehat{\tau}_{k}-\tau_{k}) =\frac{b_{k,1}}{\widehat{b}_{k,1}}\widehat{a}_{k,1}-\frac{b_{k,0}}{\widehat{b}_{k,0}}\widehat{a}_{k,0},$$

where

$$\widehat{a}_{k,1} =\frac{1}{\sqrt{n_{k}p_{k}q_{k}}\mu_{k}}\sum_{i=1}^{n_{k}}(R_{k,i}W_{k,i}-p_{k}q_{k}\mu_{k})u_{k,i}(1),$$

$$\widehat{a}_{k,0} =\frac{1}{\sqrt{n_{k}p_{k}q_{k}}(1-\mu_{k})}\sum_{i=1}^{n_{k}}(R_{k,i}(1-W_{k,i})-p_{k}q_{k}(1-\mu_{k}))u_{k,i}(0),$$

$\widehat{b}_{k,1}=(N_{k,1}\vee 1)/n_{k}$, $\widehat{b}_{k,0}=(N_{k,0}\vee 1)/n_{k}$, $b_{k,1}=p_{k}q_{k}\mu_{k}$ and $b_{k,0}=p_{k}q_{k}(1-\mu_{k})$. We will first derive the large sample distribution of

$$\widehat{a}_{k} =\widehat{a}_{k,1}-\widehat{a}_{k,0}$$

$$=\sum_{m=1}^{m_{k}}\big{(}\xi_{k,m,1}-\xi_{k,m,0}\big{)},$$

where

$$\xi_{k,m,1}=\frac{1}{\sqrt{n_{k}p_{k}q_{k}}\mu_{k}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}\big{(}R_{k,i}W_{k,i}-p_{k}q_{k}\mu_{k}\big{)}u_{k,i}(1),$$

and

$$\xi_{k,m,0}=\frac{1}{\sqrt{n_{k}p_{k}q_{k}}(1-\mu_{k})}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}\big{(}R_{k,i}(1-W_{k,i})-p_{k}q_{k}(1-\mu_{k})\big{)}u_{k,i}(0).$$

Notice that $E[\xi_{k,m,1}]=E[\xi_{k,m,0}]=0$. Moreover, notice that the terms $\xi_{k,m,1}-\xi_{k,m,0}$ are independent across clusters, $m$. In addition,

$$E[\xi_{k,m,1}^{2}]=\frac{1}{n_{k}}\sum_{i=1}^{n_{k}} 1\{m_{k,i}=m\}\frac{1-p_{k}q_{k}\mu_{k}}{\mu_{k}}u_{k,i}^{2}(1)$$

$$+\frac{2}{n_{k}}\sum_{i=1}^{n_{k}-1}\sum_{j=i+1}^{n_{k}}1\{m_{k,i}=m_{k,j}=m\}\frac{p_{k}(\sigma_{k}^{2}+\mu_{k}^{2}(1-q_{k}))}{\mu_{k}^{2}}u_{k,i}(1)u_{k,j}(1).$$

$$E[\xi_{k,m,0}^{2}]=\frac{1}{n_{k}}\sum_{i=1}^{n_{k}} 1\{m_{k,i}=m\}\frac{1-p_{k}q_{k}(1-\mu_{k})}{(1-\mu_{k})}u_{k,i}^{2}(0)$$

$$+\frac{2}{n_{k}}\sum_{i=1}^{n_{k}-1}\sum_{j=i+1}^{n_{k}}1\{m_{k,i}=m_{k,j}=m\}\frac{p_{k}(\sigma_{k}^{2}+(1-\mu_{k})^{2}(1-q_{k}))}{(1-\mu_{k})^{2}}u_{k,i}(0)u_{k,j}(0),$$

and

$$E[ \xi_{k,m,1}\xi_{k,m,0}]=-\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}p_{k}q_{k}u_{k,i}(1)u_{k,i}(0)$$

$$+\frac{1}{n_{k}}\sum_{i=1}^{n_{k}-1}\sum_{j=i+1}^{n_{k}}1\{m_{k,i}=m_{k,j}=m\}\frac{p_{k}(\mu_{k}(1-\mu_{k})(1-q_{k})-\sigma_{k}^{2})}{\mu_{k}(1-\mu_{k})}\big{(}u_{k,i}(0)u_{k,j}(1)+u_{k,i}(1)u_{k,j}(0)\big{)}.$$

We obtain:

$$n_{k}E[(\xi_{k,m,1}-\xi_{k,m,0})^{2}]$$

$$=\frac{1}{\mu_{k}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}u^{2}_{k,i}(1)+\frac{1}{1-\mu_{k}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}u^{2}_{k,i}(0)$$

$$+2p_{k}\sum_{i=1}^{n_{k}-1}\sum_{j=i+1}^{n_{k}}1\{m_{k,i}=m_{k,j}=m\}\Big{(}u_{k,i}(1)u_{k,j}(1)+u_{k,i}(0)u_{k,j}(0)-u_{k,i}(0)u_{k,j}(1)-u_{k,i}(1)u_{k,j}(0)\Big{)}$$

$$-p_{k}q_{k}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}\big{(}u^{2}_{k,i}(1)+u^{2}_{k,i}(0)-2u_{k,i}(1)u_{k,i}(0)\big{)}$$

$$\hskip 28.45274pt+2\sum_{i=1}^{n_{k}-1}\sum_{j=i+1}^{n_{k}}1\{m_{k,i}=m_{k,j}=m\}\big{(}u_{k,i}(1)u_{k,j}(1)+u_{k,i}(0)u_{k,j}(0)-u_{k,i}(0)u_{k,j}(1)-u_{k,i}(1)u_{k,j}(0)\big{)}\Bigg{)}$$

$$+2p_{k}\sigma^{2}_{k}\Bigg{(}\sum_{i=1}^{n_{k}-1}\sum_{j=i+1}^{n_{k}}1\{m_{k,i}=m_{k,j}=m\}\bigg{(}\frac{u_{k,i}(1)u_{k,j}(1)}{\mu^{2}_{k}}+\frac{u_{k,i}(0)u_{k,j}(0)}{(1-\mu_{k})^{2}}+\frac{u_{k,i}(0)u_{k,j}(1)}{\mu_{k}(1-\mu_{k})}+\frac{u_{k,i}(1)u_{k,j}(0)}{\mu_{k}(1-\mu_{k})}\Bigg{)}.$$

Therefore,

$$n_{k} E[(\xi_{k,m,1}-\xi_{k,m,0})^{2}]$$

$$=\frac{1}{\mu_{k}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}u^{2}_{k,i}(1)+\frac{1}{1-\mu_{k}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}u^{2}_{k,i}(0)$$

$$+p_{k}\Bigg{[}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}\big{(}u_{k,i}(1)-u_{k,i}(0)\big{)}\Bigg{)}^{2}-\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}\big{(}u_{k,i}(1)-u_{k,i}(0)\big{)}^{2}\Bigg{]}$$

$$-p_{k}q_{k}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}\big{(}u_{k,i}(1)-u_{k,i}(0)\big{)}\Bigg{)}^{2}$$

$$+p_{k}\sigma^{2}_{k}\Bigg{[}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}\bigg{(}\frac{u_{k,i}(1)}{\mu_{k}}+\frac{u_{k,i}(0)}{1-\mu_{k}}\bigg{)}\Bigg{)}^{2}-\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}\bigg{(}\frac{u_{k,i}(1)}{\mu_{k}}+\frac{u_{k,i}(0)}{1-\mu_{k}}\bigg{)}^{2}\Bigg{]}.$$

Let $v_{k}=\sum_{m=1}^{m_{k}}E[(\xi_{k,m,1}-\xi_{k,m,0})^{2}]$, then

$$n_{k}v_{k} =\sum_{i=1}^{n_{k}}\bigg{(}\frac{u^{2}_{k,i}(1)}{\mu_{k}}+\frac{u^{2}_{k,i}(0)}{1-\mu_{k}}\bigg{)}$$

$$+p_{k}\sum_{m=1}^{m_{k}}\Bigg{[}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}\big{(}u_{k,i}(1)-u_{k,i}(0)\big{)}\Bigg{)}^{2}-\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}\big{(}u_{k,i}(1)-u_{k,i}(0)\big{)}^{2}\Bigg{]}$$

$$-p_{k}q_{k}\sum_{m=1}^{m_{k}}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}\big{(}u_{k,i}(1)-u_{k,i}(0)\big{)}\Bigg{)}^{2}$$

$$+p_{k}\sigma^{2}_{k}\sum_{m=1}^{m_{k}}\Bigg{[}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}\bigg{(}\frac{u_{k,i}(1)}{\mu_{k}}+\frac{u_{k,i}(0)}{1-\mu_{k}}\bigg{)}\Bigg{)}^{2}-\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}\bigg{(}\frac{u_{k,i}(1)}{\mu_{k}}+\frac{u_{k,i}(0)}{1-\mu_{k}}\bigg{)}^{2}\Bigg{]}.$$

Alternatively, we can write this expression as

$$n_{k}v_{k} =\sum_{i=1}^{n_{k}}\bigg{(}\frac{u^{2}_{k,i}(1)}{\mu_{k}}+\frac{u^{2}_{k,i}(0)}{1-\mu_{k}}\bigg{)}$$

$$-p_{k}\sum_{i=1}^{n_{k}}\big{(}u_{k,i}(1)-u_{k,i}(0)\big{)}^{2}-p_{k}\sigma_{k}^{2}\sum_{i=1}^{n_{k}}\bigg{(}\frac{u_{k,i}(1)}{\mu_{k}}+\frac{u_{k,i}(0)}{1-\mu_{k}}\bigg{)}^{2}$$

$$+p_{k}(1-q_{k})\sum_{m=1}^{m_{k}}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}\big{(}u_{k,i}(1)-u_{k,i}(0)\big{)}\Bigg{)}^{2}$$

$$+p_{k}\sigma^{2}_{k}\sum_{m=1}^{m_{k}}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}\bigg{(}\frac{u_{k,i}(1)}{\mu_{k}}+\frac{u_{k,i}(0)}{1-\mu_{k}}\bigg{)}\Bigg{)}^{2}.$$

The sum of the first three terms is minimized for $p_{k}=1$ and $\sigma_{k}^{2}=\mu_{k}(1-\mu_{k})$, in which case this sum is equal to zero. Therefore,

$$v_{k} \geq(p_{k}\min_{m}n_{k,m})(1-q_{k})\sum_{m=1}^{m_{k}}\frac{n_{k,m}}{n_{k}}\Bigg{(}\frac{1}{n_{k,m}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}\big{(}u_{k,i}(1)-u_{k,i}(0)\big{)}\Bigg{)}^{2}$$

$$+(p_{k}\min_{m}n_{k,m})\,\sigma^{2}_{k}\sum_{m=1}^{m_{k}}\frac{n_{k,m}}{n_{k}}\Bigg{(}\frac{1}{n_{k,m}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}\bigg{(}\frac{u_{k,i}(1)}{\mu_{k}}+\frac{u_{k,i}(0)}{1-\mu_{k}}\bigg{)}\Bigg{)}^{2}.$$ (A.1)

We will assume that $\liminf_{k\rightarrow\infty}((1-q_{k})\vee\sigma_{k}^{2})>0$, so either sampling or assignment or both are correlated within cluster. (We study the case $q_{k}=1$ and $\sigma_{k}^{2}=0$ separately below.) In addition, assume (i) $\liminf_{k\rightarrow\infty}(1-q_{k})>0$ and

$$\liminf\limits_{k\rightarrow\infty}\sum_{m=1}^{m_{k}}\frac{n_{k,m}}{n_{k}}\Bigg{(}\frac{1}{n_{k,m}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}\big{(}u_{k,i}(1)-u_{k,i}(0)\big{)}\Bigg{)}^{2}>0,$$ (A.2)

or (ii) $\liminf_{k\rightarrow\infty}\sigma_{k}^{2}>0$ and

$$\liminf\limits_{k\rightarrow\infty}\sum_{m=1}^{m_{k}}\frac{n_{k,m}}{n_{k}}\Bigg{(}\frac{1}{n_{k,m}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}\bigg{(}\frac{u_{k,i}(1)}{\mu_{k}}+\frac{u_{k,i}(0)}{1-\mu_{k}}\bigg{)}\Bigg{)}^{2}>0.$$ (A.3)

Equation (A.2) would be violated if, as $k$ increases, there is no variation in average treatment effects across clusters. Equation (A.3) would be violated if as $k$ increases there is no variation in average potential outcomes across clusters. If equations (A.2) and (A.3) hold, $v_{k}$ is bounded below by a term of order at least $p_{k}\min_{m}n_{k,m}$. Recall our assumption, $\liminf_{k\rightarrow\infty}p_{k}\min_{m}n_{k,m}>0$, so the average number of observations sampled per cluster, conditional on the cluster being sampled, does not go to zero. Then,

$$\liminf\limits_{k\rightarrow\infty}v_{k}>0.$$

To obtain a CLT, we will check Lyapunov’s condition,

$$\lim_{k\rightarrow\infty}\sum_{m=1}^{m_{k}}\frac{1}{v_{k}^{1+\delta/2}}E[|\xi_{k,m,1}-\xi_{k,m,0}|^{2+\delta}]=0,$$

for some $\delta>0$. Because potential outcomes are uniformly bounded and $\mu_{k}$ is uniformly bounded away from zero, we obtain

$$|\xi_{k,m,1}|^{2+\delta} \leq c\frac{n_{k,m}^{2+\delta}}{(n_{k}p_{k}q_{k})^{1+\delta/2}}\left|\frac{1}{n_{k,m}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}|R_{k,i}W_{k,i}-p_{k}q_{k}\mu_{k}|\right|^{2+\delta},$$

where $c$ is some generic positive constant, whose value may change across equations. Consider $\delta=1$, and let

$$S {}_{k,m,1}^{3}=E\left[\left|\frac{1}{n_{k,m}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}|R_{k,i}W_{k,i}-p_{k}q_{k}\mu_{k}|\right|^{3}\right]$$

$$\leq\frac{1}{n_{k,m}^{3}}n_{k,m}E[|R_{k,i}W_{k,i}-p_{k}q_{k}\mu_{k}|^{3}]$$

$$+\frac{3}{n_{k,m}^{3}}n_{k,m}(n_{k,m}-1)E\big{[}|R_{k,i}W_{k,i}-p_{k}q_{k}\mu_{k}|^{2}|R_{k,j}W_{k,j}-p_{k}q_{k}\mu_{k}|\big{|}m_{k,i}=m_{k,j}=m\big{]}$$

$$+\frac{6}{n_{k,m}^{3}}\dbinom{n_{k,m}}{3}E\big{[}|R_{k,i}W_{k,i}-p_{k}q_{k}\mu_{k}||R_{k,j}W_{k,j}-p_{k}q_{k}\mu_{k}||R_{k,t}W_{k,t}-p_{k}q_{k}\mu_{k}|\big{|}m_{k,i}=m_{k,j}=m_{k,t}=m\big{]},$$

for $i\neq j\neq t$. (The second and third terms on the left-hand side of last equation only appear when $n_{k,m}\geq 2$ and $n_{k,m}\geq 3$, respectively) As a result,

$$S_{k,m,1}^{3} \leq c\left(\frac{p_{k}q_{k}}{n_{k,m}^{2}}+\frac{p_{k}^{2}q_{k}}{n_{k,m}}+p_{k}^{3}q_{k}\right)$$

$$\leq c\,p_{k}^{3}q_{k}\left(\frac{1}{p_{k}^{2}\min_{m}n_{k,m}^{2}}+\frac{1}{p_{k}\min_{m}n_{k,m}}+1\right).$$

Because $\liminf_{k\rightarrow\infty}p_{k}\min_{m}n_{k,m}>0$, for large enough $k$ we obtain,

$$E[|\xi_{k,m,1}|^{3}] \leq c\frac{p_{k}^{3}q_{k}n_{k,m}^{3}}{(n_{k}p_{k}q_{k})^{3/2}},$$

and the same bound applies for $E[|\xi_{k,m,0}|^{3}]$. Notice that

$$\sum_{m=1}^{m_{k}}E[|\xi_{k,m,1}-\xi_{k,m,0}|^{3}] \leq\sum_{m=1}^{m_{k}}E[(|\xi_{k,m,1}|+|\xi_{k,m,0}|)^{3}]$$

$$=\sum_{m=1}^{m_{k}}E[|\xi_{k,m,1}|^{3}]+\sum_{m=1}^{m_{k}}E[|\xi_{k,m,0}|^{3}]$$

$$+3\sum_{m=1}^{m_{k}}E[|\xi_{k,m,1}|^{2}|\xi_{k,m,0}|]+3\sum_{m=1}^{m_{k}}E[|\xi_{k,m,1}||\xi_{k,m,0}|^{2}].$$

Now, Hölder’s inequality implies that

$$\frac{p_{k}^{3}q_{k}\sum_{m=1}^{m_{k}}n_{k,m}^{3}}{v_{k}^{3/2}(n_{k}p_{k}q_{k})^{3/2}}\longrightarrow 0,$$ (A.4)

is sufficient for the Lyapunov condition to hold. Because $\max_{m}n_{k,m}/\min_{m}n_{k,m}$ is bounded asymptotically, we obtain,

$$\limsup\limits_{k\rightarrow\infty}\frac{p_{k}^{3}q_{k}\sum_{m=1}^{m_{k}}n_{k,m}^{3}}{v_{k}^{3/2}(n_{k}p_{k}q_{k})^{3/2}} \leq\limsup\limits_{k\rightarrow\infty}c\,\frac{p_{k}^{3}q_{k}m_{k}\max_{m}n_{k,m}^{3}}{(p_{k}^{2}q_{k}m_{k}\min_{m}n_{k,m}^{2})^{3/2}}$$

$$\leq\limsup\limits_{k\rightarrow\infty}\left(\frac{\max_{m}n_{k,m}}{\min_{m}n_{k,m}}\right)^{3}\frac{c}{\sqrt{q_{k}m_{k}}}=0,$$

and so the Lyapunov condition holds. As a result, we obtain

$$\widehat{a}_{k}/\sqrt{v_{k}}\stackrel{{\scriptstyle d}}{{\longrightarrow}}N(0,1).$$

We will next prove that both $\widehat{a}_{k,1}/\sqrt{v_{k}}$ and $\widehat{a}_{k,0}/\sqrt{v_{k}}$ are $\mathcal{O}_{p}(1)$.

$$E[\widehat{a}_{k,1}^{2}] =\frac{1}{n_{k}p_{k}q_{k}}\frac{1}{\mu_{k}^{2}}\sum_{m=1}^{m_{k}}E\left[\left(\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}(R_{k,i}W_{k,i}-p_{k}q_{k}\mu_{k})u_{k,i}(1)\right)^{2}\right]$$

$$\leq c\,\frac{1}{n_{k}p_{k}q_{k}}\sum_{m=1}^{m_{k}}\big{(}n_{k,m}p_{k}q_{k}+n_{k,m}(n_{k,m}-1)p_{k}^{2}q_{k}\big{)}$$

$$=c\,\Big{(}1+\sum_{m=1}^{m_{k}}\frac{n_{k,m}(n_{k,m}-1)p_{k}}{n_{k}}\Big{)}.$$

Therefore,

$$E[(\widehat{a}_{k,1}/\sqrt{v_{k}})^{2}] \leq c\left(\frac{1}{p_{k}\min_{m}n_{k,m}}+\sum_{m=1}^{m_{k}}\frac{(\max_{m}n_{k,m})(n_{k,m}-1)p_{k}}{n_{k}p_{k}\min_{m}n_{k,m}}\right).$$

Because $\limsup\max_{m}n_{k,m}/\min_{m}n_{k,m}<\infty$, we obtain $\limsup_{k\rightarrow\infty}E[(\widehat{a}_{k,1}/v_{k})^{2}]<\infty$. As a result, $\widehat{a}_{k,1}/\sqrt{v_{k}}$ is $\mathcal{O}_{p}(1)$.

Let $\widetilde{b}_{k,1}=N_{k,1}/n_{k}$. Consider $k$ large enough, so $p_{k}\min_{m}n_{k,m}$ is bounded away from zero, making $\widetilde{b}_{k,1}/b_{k,1}$ well-defined. Notice that $E[\widetilde{b}_{k,1}/b_{k,1}]=1$ and

$$\mbox{var}(\widetilde{b}_{k,1}/b_{k,1}) =\frac{1}{(n_{k}p_{k}q_{k}\mu_{k})^{2}}\sum_{m=1}^{m_{k}}E\left[\left(\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}(R_{k,i}W_{k,i}-n_{k}p_{k}q_{k}\mu_{k})\right)^{2}\right]$$

$$=\frac{n_{k}p_{k}q_{k}\mu_{k}(1-p_{k}q_{k}\mu_{k})}{(n_{k}p_{k}q_{k}\mu_{k})^{2}}+\sum_{m=1}^{m_{k}}\frac{n_{k,m}(n_{k,m}-1)p_{k}^{2}q_{k}(\sigma_{k}^{2}+(1-q_{k})\mu_{k}^{2})}{(n_{k}p_{k}q_{k}\mu_{k})^{2}}$$

$$\leq\frac{1-p_{k}q_{k}\mu_{k}}{n_{k}p_{k}q_{k}\mu_{k}}+c\,\frac{n_{k}(\max_{m}n_{k,m}-1)p_{k}^{2}q_{k}}{(n_{k}p_{k}q_{k})^{2}}$$

$$\leq\frac{1-p_{k}q_{k}\mu_{k}}{n_{k}p_{k}q_{k}\mu_{k}}+c\,\frac{(\max_{m}n_{k,m}-1)}{\min n_{k,m}}\frac{1}{q_{k}m_{k}}\longrightarrow 0.$$

This implies $\widetilde{b}_{k,1}/b_{k,1}\stackrel{{\scriptstyle p}}{{\rightarrow}}1$. Analogous calculations yield $\widetilde{b}_{0,k}/b_{0,k}\stackrel{{\scriptstyle p}}{{\rightarrow}}1$. For large enough $k$, $\widetilde{b}_{k,1}/b_{k,1}=0$ if and only if $N_{k,1}=0$, which implies $\Pr(N_{k,1}=0)\rightarrow 0$. It follows that, for large enough $k$,

$$\Pr(|\widetilde{b}_{k,1}/b_{k,1}-\widehat{b}_{k,1}/b_{k,1}|=0)=\Pr(N_{k,1}>0)\longrightarrow 1$$

and $\widehat{b}_{k,1}/b_{k,1}\stackrel{{\scriptstyle p}}{{\rightarrow}}1$. Using analogous calculations, we obtain $\widehat{b}_{k,0}/b_{k,0}\stackrel{{\scriptstyle p}}{{\rightarrow}}1$. As a result,

$$\sqrt{n_{k}p_{k}q_{k}}(\widehat{\tau}_{k}-\tau_{k})/v_{k}^{1/2} =\frac{b_{k,1}}{\widehat{b}_{k,1}}\frac{\widehat{a}_{k,1}}{v_{k}^{1/2}}-\frac{b_{k,0}}{\widehat{b}_{k,0}}\frac{\widehat{a}_{k,0}}{v_{k}^{1/2}}$$

$$=\frac{\widehat{a}_{k}}{v_{k}^{1/2}}+\left(\frac{b_{k,1}}{\widehat{b}_{k,1}}-1\right)\frac{\widehat{a}_{k,1}}{v_{k}^{1/2}}-\left(\frac{b_{k,0}}{\widehat{b}_{k,0}}-1\right)\frac{\widehat{a}_{k,0}}{v_{k}^{1/2}}$$

$$=\widehat{a}_{k}/\sqrt{v_{k}}+\scaleto{\mathcal{O}}{5pt}_{p}(1).$$

Therefore,

$$\sqrt{n_{k}p_{k}q_{k}}(\widehat{\tau}_{k}-\tau_{k})/v_{k}^{1/2}\stackrel{{\scriptstyle d}}{{\longrightarrow}}N(0,1).$$

Using $\widetilde{b}_{1,k}/b_{1,k}\stackrel{{\scriptstyle p}}{{\rightarrow}}1$ and $\widetilde{b}_{0,k}/b_{0,k}\stackrel{{\scriptstyle p}}{{\rightarrow}}1$, it is easy to show $N_{k}/(n_{k}p_{k}q_{k})\stackrel{{\scriptstyle p}}{{\rightarrow}}1$, which implies

$$\sqrt{N_{k}}(\widehat{\tau}_{k}-\tau_{k})/v_{k}^{1/2}\stackrel{{\scriptstyle d}}{{\longrightarrow}}N(0,1).$$

We will next consider the case of $q_{k}=1$ and $\sigma_{k}^{2}=0$, where no clustering is required. Consider

$$\vartheta_{k,i,1}=\frac{1}{\sqrt{n_{k}p_{k}}\mu_{k}}\big{(}R_{k,i}W_{k,i}-p_{k}\mu_{k}\big{)}u_{k,i}(1)$$

and

$$\vartheta_{k,i,0}=\frac{1}{\sqrt{n_{k}p_{k}}(1-\mu_{k})}\big{(}R_{k,i}(1-W_{k,i})-p_{k}(1-\mu_{k})\big{)}u_{k,i}(0).$$

Redefine now $v_{k}=\sum_{i=1}^{n_{k}}E\big{[}(\vartheta_{k,i,1}-\vartheta_{k,i,0})^{2}\big{]}$. Then,

$$v_{k} =\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}\bigg{(}\frac{u^{2}_{k,i}(1)}{\mu_{k}}+\frac{u^{2}_{k,i}(0)}{1-\mu_{k}}\bigg{)}-p_{k}\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}\big{(}u_{k,i}(1)-u_{k,i}(0)\big{)}^{2}.$$

Notice that $v_{k}$ is minimized for $p_{k}=1$, in which case

$$v_{k} =\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}\bigg{(}\frac{u^{2}_{k,i}(1)}{\mu_{k}}+\frac{u^{2}_{k,i}(0)}{1-\mu_{k}}\bigg{)}-\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}\big{(}u_{k,i}(1)-u_{k,i}(0)\big{)}^{2}$$

$$=\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}\bigg{(}\frac{1-\mu_{k}}{\mu_{k}}u^{2}_{k,i}(1)+\frac{\mu_{k}}{1-\mu_{k}}u^{2}_{k,i}(0)+2u_{k,i}(1)u_{k,i}(0)\bigg{)}$$

$$=\mu_{k}(1-\mu_{k})\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}\bigg{(}\frac{u^{2}_{k,i}(1)}{\mu^{2}_{k}}+\frac{u^{2}_{k,i}(0)}{(1-\mu_{k})^{2}}+2\,\frac{u_{k,i}(1)u_{k,i}(0)}{\mu_{k}(1-\mu_{k})}\bigg{)}$$

$$=\mu_{k}(1-\mu_{k})\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}\bigg{(}\frac{u_{k,i}(1)}{\mu_{k}}+\frac{u_{k,i}(0)}{1-\mu_{k}}\bigg{)}^{2}.$$

Therefore, the assumption

$$\liminf\limits_{k\rightarrow\infty}\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}\bigg{(}\frac{u_{k,i}(1)}{\mu_{k}}+\frac{u_{k,i}(0)}{1-\mu_{k}}\bigg{)}^{2}>0$$

is enough for $\liminf_{k\rightarrow\infty}v_{k}>0$. Notice now that

$$E[|\vartheta_{k,i,1}|^{3}] \leq\frac{1}{(n_{k}p_{k})^{3/2}}E[|R_{k,i}W_{k,i}-p_{k}\mu_{k}|^{3}]$$

$$=\frac{1}{(n_{k}p_{k})^{3/2}}(1-p_{k}\mu_{k})^{3}p_{k}\mu_{k}+(p_{k}\mu_{k})^{3}(1-p_{k}\mu_{k})$$

$$\leq c\frac{p_{k}}{(n_{k}p_{k})^{3/2}},$$

and the same bound holds for $E[|\vartheta_{k,i,0}|^{3}]$. Therefore, for the Lyapunov condition to hold, it is enough that

$$\frac{n_{k}p_{k}}{(n_{k}p_{k})^{3/2}}=\frac{1}{\sqrt{n_{k}p_{k}}}\longrightarrow 0,$$

or $n_{k}p_{k}\rightarrow\infty$. That is, assumptions (i)-(iii), which we used for the clustered case, are replaced by $n_{k}p_{k}\rightarrow\infty$.

A.2.2 Estimation of the variance

Let $\widehat{U}_{k,i}=Y_{k,i}-\widehat{\alpha}_{k}-\widehat{\tau}_{k}W_{k,i}$ be the residuals from the regression of $Y_{k,i}$ or a constant and $W_{k,i}$. Here, $\widehat{\alpha}_{k}$ is the coefficient on the constant regressor equal to one, and $\widehat{\tau}_{k}$ is the coefficient on $W_{k,i}$. We have already shown $v_{k}^{-1/2}(\widehat{\tau}_{k}-\tau_{k})=\mathcal{O}_{p}(1/\sqrt{n_{k}p_{k}q_{k}})$. The same is true about $\widehat{\alpha}_{k}$ (e.g., apply the proof for $\widehat{\tau}_{k}$ after replacing each $y_{k,i}(1)$ with a zero). Define $\widehat{\Sigma}_{k}=\sum_{m=1}^{m_{k}}\widehat{\Sigma}_{k,m}$, where

$$\widehat{\Sigma}_{k,m}=\left(\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}\left(\begin{array}[]{c}\widehat{U}_{k,i}\\ W_{k,i}\widehat{U}_{k,i}\end{array}\right)\right)\left(\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}\left(\begin{array}[]{c}\widehat{U}_{k,i}\\ W_{k,i}\widehat{U}_{k,i}\end{array}\right)\right)^{\prime}.$$

Also, let

$$\widehat{Q}_{k}=\sum_{i=1}^{n_{k}}R_{k,i}\left(\begin{array}[]{c}1\\ W_{k,i}\end{array}\right)\left(\begin{array}[]{c}1\\ W_{k,i}\end{array}\right)^{\prime},$$

and $z=(0,1)^{\prime}$. Then, the cluster estimator of the variance of $\sqrt{N_{k}}(\widehat{\tau}_{k}-\tau_{k})$ is

$$\widehat{V}_{k}^{\rm{cluster}}=N_{k}z^{\prime}\widehat{Q}_{k}^{-1}\widehat{\Sigma}_{k}\widehat{Q}_{k}^{-1}z.$$

Notice that

$$(n_{k}p_{k}q_{k})^{-1}E[\widehat{Q}_{k}]=\left(\begin{array}[]{cc}1&\mu_{k}\\ \mu_{k}&\mu_{k}\end{array}\right).$$

In addition,

$$\frac{1}{n_{k}p_{k}q_{k}}\widehat{Q}_{k}(2,2)=\frac{1}{n_{k}p_{k}q_{k}}\sum_{m=1}^{m_{k}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}.$$

$$\mbox{var}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}\Bigg{)} =n_{k,m}p_{k}q_{k}\mu_{k}(1-p_{k}q_{k}\mu_{k})$$

$$+n_{k,m}(n_{k,m}-1)p_{k}^{2}q_{k}(\sigma_{k}^{2}+\mu_{k}^{2}(1-q_{k})).$$

Therefore, under conditions (i)-(iii), we obtain

$$\mbox{var}\Bigg{(}\frac{1}{n_{k}p_{k}q_{k}}\widehat{Q}_{k}(2,2)\Bigg{)} \leq\frac{c}{n_{k}p_{k}q_{k}}\Bigg{(}1+p_{k}(\textstyle\max_{m}n_{k,m}-1)\Bigg{)}$$

$$=c\,\frac{\textstyle\max_{m}n_{k,m}}{n_{k}q_{k}}+\scaleto{\mathcal{O}}{5pt}(1)$$

$$\leq c\,\frac{\textstyle\max_{m}n_{k,m}}{\textstyle\min_{m}n_{k,m}}\,\frac{1}{q_{k}m_{k}}+\scaleto{\mathcal{O}}{5pt}(1)\longrightarrow 0.$$

Analogous calculations yield $\mbox{var}((n_{k}p_{k}q_{k})^{-1}\widehat{Q}_{k}(1,1))\rightarrow 0$. Therefore,

$$\frac{1}{n_{k}p_{k}q_{k}}\widehat{Q}_{k}=\left(\begin{array}[]{cc}1&\mu_{k}\\ \mu_{k}&\mu_{k}\end{array}\right)+\scaleto{\mathcal{O}}{5pt}_{p}(1)$$

and

$$n_{k}q_{k}p_{k}\widehat{Q}_{k}^{-1}=H_{k}+\scaleto{\mathcal{O}}{5pt}_{p}(1),\quad\mbox{where}\quad H_{k}=\frac{1}{\mu_{k}(1-\mu_{k})}\left(\begin{array}[]{rc}\mu_{k}&-\mu_{k}\\ -\mu_{k}&\ \ 1\end{array}\right).$$

Now, let $U_{k,i}=Y_{k,i}-\alpha_{k}-\tau_{k}W_{k,i}=W_{k,i}u_{k,i}(1)+(1-W_{k,i})u_{k,i}(0)$. Notice that

$$v_{k}^{-1/2}\max_{i=1,\ldots,n_{k}}|\widehat{U}_{k,i}-U_{k,i}|\leq v_{k}^{-1/2}|\widehat{\alpha}_{k}-\alpha_{k}|+v_{k}^{-1/2}|\widehat{\tau}_{k}-\tau_{k}|=\mathcal{O}_{p}(1/\sqrt{n_{k}p_{k}q_{k}}).$$

Define $\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\Sigma$\kern-0.18004pt}}}_{k}=\sum_{m=1}^{m_{k}}\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\Sigma$\kern-0.18004pt}}}_{k,m}$, where

$$\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\Sigma$\kern-0.18004pt}}}_{k,m}=\left(\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}\left(\begin{array}[]{c}U_{k,i}\\ W_{k,i}U_{k,i}\end{array}\right)\right)\left(\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}\left(\begin{array}[]{c}U_{k,i}\\ W_{k,i}U_{k,i}\end{array}\right)\right)^{\prime}.$$

We will show

$$\frac{1}{n_{k}p_{k}q_{k}v_{k}}(\widehat{\Sigma}_{k}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\Sigma$\kern-0.18004pt}}}_{k})\stackrel{{\scriptstyle p}}{{\longrightarrow}}0.$$

Notice that

$$\widehat{\Sigma}_{k,m}(2,2)-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\Sigma$\kern-0.18004pt}}}_{k,m}(2,2) =\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}(\widehat{U}_{k,i}-U_{k,i})\Bigg{)}^{2}$$

$$+2\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}U_{k,i}\Bigg{)}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}(\widehat{U}_{k,i}-U_{k,i})\Bigg{)}.$$

Therefore,

$$\frac{1}{n_{k}p_{k}q_{k}v_{k}}|\widehat{\Sigma}_{k}(2,2)-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\Sigma$\kern-0.18004pt}}}_{k}(2,2)| \leq c\frac{1}{n_{k}p_{k}q_{k}v_{k}}\sum_{m=1}^{m_{k}}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}\Bigg{)}^{2}$$

$$\qquad\times\Bigg{(}\max_{i=1,\ldots,n_{k}}|\widehat{U}_{k,i}-U_{k,i}|^{2}+\max_{i=1,\ldots,n_{k}}|\widehat{U}_{k,i}-U_{k,i}|\Bigg{)}.$$

The same expression holds for the off-diagonal elements of $\widehat{\Sigma}_{k,m}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\Sigma$\kern-0.18004pt}}}_{k,m}$. For $\widehat{\Sigma}_{k,m}(1,1)-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\Sigma$\kern-0.18004pt}}}_{k,m}(1,1)$, the expression holds once we replace each $W_{k,i}$ with a one. Let $\|\cdot\|$ be the Frobenius norm of a matrix. Then,

$$\frac{1}{n_{k}p_{k}q_{k}v_{k}}\|\widehat{\Sigma}_{k}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\Sigma$\kern-0.18004pt}}}_{k}\| \leq c\frac{1}{n_{k}p_{k}q_{k}v_{k}}\sum_{m=1}^{m_{k}}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}\Bigg{)}^{2}$$

$$\qquad\times\Bigg{(}\max_{i=1,\ldots,n_{k}}|\widehat{U}_{k,i}-U_{k,i}|^{2}+\max_{i=1,\ldots,n_{k}}|\widehat{U}_{k,i}-U_{k,i}|\Bigg{)}.$$

We will prove that the right-hand side of the previous equation converges to zero in probability. We will factorize each term into a expression that is bounded in probability and one that converges to zero in $L_{1}$.

$$E\Bigg{[}\sum_{m=1}^{m_{k}}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}\Bigg{)}^{2}\Bigg{]} \leq n_{k}p_{k}q_{k}+n_{k}(\max_{m}n_{k,m}-1)p_{k}^{2}q_{k}.$$

For the first term, notice that

$$\max_{i=1,\ldots,n_{k}}|\widehat{U}_{k,i} -U_{k,i}|^{2}\frac{n_{k}p_{k}q_{k}+n_{k}(\max_{m}n_{k,m}-1)p_{k}^{2}q_{k}}{n_{k}p_{k}q_{k}v_{k}}$$

$$=\frac{n_{k}p_{k}q_{k}}{v_{k}}\max_{i=1,\ldots,n_{k}}|\widehat{U}_{k,i}-U_{k,i}|^{2}\Bigg{(}\frac{n_{k}p_{k}q_{k}+n_{k}(\max_{m}n_{k,m}-1)p_{k}^{2}q_{k}}{(n_{k}p_{k}q_{k})^{2}}\Bigg{)}$$

$$\leq\frac{n_{k}p_{k}q_{k}}{v_{k}}\max_{i=1,\ldots,n_{k}}|\widehat{U}_{k,i}-U_{k,i}|^{2}\Bigg{(}\frac{1}{n_{k}p_{k}q_{k}}+\frac{\max_{m}n_{k,m}-1}{\min_{m}n_{k,m}}\frac{1}{q_{k}m_{k}}\Bigg{)}$$

$$=\mathcal{O}_{p}(1)\,\scaleto{\mathcal{O}}{5pt}(1).$$

For the second term, using the fact that $v_{k}$ is greater or equal to $p_{k}\min_{m}n_{k,m}>0$ times a term with limit inferior that is bounded away from zero, we obtain

$$\max_{i=1,\ldots,n_{k}}|\widehat{U}_{k,i} -U_{k,i}|\frac{n_{k}p_{k}q_{k}+n_{k}(\max_{m}n_{k,m}-1)p_{k}^{2}q_{k}}{n_{k}p_{k}q_{k}v_{k}}$$

$$=\Big{(}\frac{n_{k}p_{k}q_{k}}{v_{k}}\Big{)}^{1/2}\!\!\max_{i=1,\ldots,n_{k}}|\widehat{U}_{k,i}-U_{k,i}|\Bigg{(}\frac{n_{k}p_{k}q_{k}+n_{k}(\max_{m}n_{k,m}-1)p_{k}^{2}q_{k}}{(n_{k}p_{k}q_{k})^{3/2}v_{k}^{1/2}}\Bigg{)}$$

$$\leq\Big{(}\frac{n_{k}p_{k}q_{k}}{v_{k}}\Big{)}^{1/2}\!\!\max_{i=1,\ldots,n_{k}}|\widehat{U}_{k,i}-U_{k,i}|\Bigg{(}\frac{1}{(n_{k}p_{k}q_{k}v_{k})^{1/2}}+\frac{\max_{m}n_{k,m}-1}{\min_{m}n_{k,m}}\frac{1}{(q_{k}m_{k})^{1/2}}\Bigg{)}$$

$$=\mathcal{O}_{p}(1)\,\scaleto{\mathcal{O}}{5pt}(1).$$

As a result, we obtain

$$\frac{1}{n_{k}p_{k}q_{k}v_{k}}\|\widehat{\Sigma}_{k}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\Sigma$\kern-0.18004pt}}}_{k}\|=\scaleto{\mathcal{O}}{5pt}_{p}(1).$$

Notice that

$$\frac{n_{k}p_{k}q_{k}}{v_{k}} \widehat{Q}_{k}^{-1}\widehat{\Sigma}_{k}\widehat{Q}_{k}^{-1}-H_{k}\frac{\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\Sigma$\kern-0.18004pt}}}_{k}}{n_{k}p_{k}q_{k}v_{k}}H_{k}=H_{k}\frac{\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\Sigma$\kern-0.18004pt}}}_{k}}{n_{k}p_{k}q_{k}v_{k}}\Big{(}n_{k}p_{k}q_{k}\widehat{Q}_{k}^{-1}-H_{k}\Big{)}$$

$$+\Big{(}n_{k}p_{k}q_{k}\widehat{Q}_{k}^{-1}-H_{k}\Big{)}\frac{\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\Sigma$\kern-0.18004pt}}}_{k}}{n_{k}p_{k}q_{k}v_{k}}\Big{(}n_{k}p_{k}q_{k}\widehat{Q}_{k}^{-1}\Big{)}+\Big{(}n_{k}p_{k}q_{k}\widehat{Q}_{k}^{-1}\Big{)}\frac{\widehat{\Sigma}_{k}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\Sigma$\kern-0.18004pt}}}_{k}}{n_{k}p_{k}q_{k}v_{k}}\Big{(}n_{k}p_{k}q_{k}\widehat{Q}_{k}^{-1}\Big{)}.$$

Therefore, to show that the left-hand side of the last equation is $\scaleto{\mathcal{O}}{5pt}_{p}(1)$, it is only left to show that $\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\Sigma$\kern-0.18004pt}}}_{k}/(n_{k}p_{k}q_{k}v_{k})$ is $\mathcal{O}_{p}(1)$. We will prove this next. Notice that

$$\frac{1}{n_{k}p_{k}q_{k}v_{k}}\|\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\Sigma$\kern-0.18004pt}}}_{k}\| \leq c\frac{1}{n_{k}p_{k}q_{k}v_{k}}\sum_{m=1}^{m_{k}}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}\Bigg{)}^{2}.$$

Therefore,

$$E\Bigg{[}\frac{1}{n_{k}p_{k}q_{k}v_{k}}\|\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\Sigma$\kern-0.18004pt}}}_{k}\|\Bigg{]} \leq c\frac{1}{n_{k}p_{k}q_{k}v_{k}}\Bigg{(}n_{k}p_{k}q_{k}+n_{k}(\max_{m}n_{k,m}-1)p_{k}^{2}q_{k}\Bigg{)}.$$

Then,

$$E\Bigg{[}\frac{1}{n_{k}p_{k}q_{k}v_{k}}\|\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\Sigma$\kern-0.18004pt}}}_{k}\|\Bigg{]} \leq c\Bigg{(}\frac{1}{v_{k}}+\frac{p_{k}(\max_{m}n_{k,m}-1)}{p_{k}\min_{m}n_{k,m}}\Bigg{)}<\infty.$$

We, therefore, obtain,

$$\frac{n_{k}p_{k}q_{k}}{v_{k}}\widehat{Q}_{k}^{-1}\widehat{\Sigma}_{k}\widehat{Q}_{k}^{-1}-H_{k}\frac{\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\Sigma$\kern-0.18004pt}}}_{k}}{n_{k}p_{k}q_{k}v_{k}}H_{k}\stackrel{{\scriptstyle p}}{{\longrightarrow}}0.$$

Because $N_{k}/(n_{k}p_{k}q_{k})\stackrel{{\scriptstyle p}}{{\rightarrow}}1$, we obtain

$$\widehat{V}_{k}^{\rm{cluster}}/v_{k} =z^{\prime}H_{k}\frac{\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\Sigma$\kern-0.18004pt}}}_{k}}{n_{k}p_{k}q_{k}v_{k}}H_{k}z+\scaleto{\mathcal{O}}{5pt}_{p}(1)$$

$$=\frac{1}{n_{k}p_{k}q_{k}v_{k}}\left(\frac{1}{\mu_{k}(1-\mu_{k})}\right)^{2}\sum_{m=1}^{m_{k}}\bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}(W_{k,i}-\mu_{k})U_{k,i}\bigg{)}^{2}+\scaleto{\mathcal{O}}{5pt}_{p}(1).$$

Recall that $U^{2}_{k,i}=u^{2}_{k,i}(1)W_{k,i}+u^{2}_{k,i}(0)(1-W_{k,i})$. Notice that

$$E\bigg{[}\bigg{(}\sum_{i=1}^{n_{k}} 1\{m_{k,i}=m\}R_{k,i}(W_{k,i}-\mu_{k})U_{k,i}\bigg{)}^{2}\bigg{]}$$

$$=\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}p_{k}q_{k}\mu_{k}(1-\mu_{k})\Big{(}(1-\mu_{k})u^{2}_{k,i}(1)+\mu_{k}u^{2}_{k,i}(0)\Big{)}$$

$$+2\sum_{i=1}^{n_{k}-1}\sum_{j=i+1}^{n_{k}}1\{m_{k,i}=m_{k,j}=m\}p_{k}^{2}q_{k}\Big{[}(\sigma_{k}^{2}+\mu_{k}^{2})(1-\mu_{k})^{2}u_{k,i}(1)u_{k,j}(1)$$

$$\qquad+\mu_{k}(1-\mu_{k})(\sigma_{k}^{2}-\mu_{k}(1-\mu_{k}))(u_{k,i}(0)u_{k,j}(1)+u_{k,i}(1)u_{k,j}(0))$$

$$\qquad+(\sigma_{k}^{2}+(1-\mu_{k})^{2})\mu_{k}^{2}u_{k,i}(0)u_{k,j}(0)\Big{]}.$$

Let

$$v_{k}^{\rm cluster}=\frac{1}{n_{k}p_{k}q_{k}}\left(\frac{1}{\mu_{k}(1-\mu_{k})}\right)^{2}\sum_{m=1}^{m_{k}}E\bigg{[}\bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}(W_{k,i}-\mu_{k})U_{k,i}\bigg{)}^{2}\bigg{]}.$$

Then,

$$n_{k}v_{k}^{\rm{cluster}} =\sum_{i=1}^{n_{k}}\bigg{(}\frac{u^{2}_{k,i}(1)}{\mu_{k}}+\frac{u^{2}_{k,i}(0)}{1-\mu_{k}}\bigg{)}$$

$$+p_{k}\sum_{m=1}^{m_{k}}\Bigg{[}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}\big{(}u_{k,i}(1)-u_{k,i}(0)\big{)}\Bigg{)}^{2}-\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}\big{(}u_{k,i}(1)-u_{k,i}(0)\big{)}^{2}\Bigg{]}$$

$$+p_{k}\sigma^{2}_{k}\sum_{m=1}^{m_{k}}\Bigg{[}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}\bigg{(}\frac{u_{k,i}(1)}{\mu_{k}}+\frac{u_{k,i}(0)}{1-\mu_{k}}\bigg{)}\Bigg{)}^{2}-\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}\bigg{(}\frac{u_{k,i}(1)}{\mu_{k}}+\frac{u_{k,i}(0)}{1-\mu_{k}}\bigg{)}^{2}\Bigg{]}.$$

Alternatively, we can write

$$n_{k}v_{k}^{\rm cluster} =\sum_{i=1}^{n_{k}}\bigg{(}\frac{u^{2}_{k,i}(1)}{\mu_{k}}+\frac{u^{2}_{k,i}(0)}{1-\mu_{k}}\bigg{)}$$

$$-p_{k}\sum_{i=1}^{n_{k}}\big{(}u_{k,i}(1)-u_{k,i}(0)\big{)}^{2}-p_{k}\sigma_{k}^{2}\sum_{i=1}^{n_{k}}\bigg{(}\frac{u_{k,i}(1)}{\mu_{k}}+\frac{u_{k,i}(0)}{1-\mu_{k}}\bigg{)}^{2}$$

$$+p_{k}\sum_{m=1}^{m_{k}}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}\big{(}u_{k,i}(1)-u_{k,i}(0)\big{)}\Bigg{)}^{2}$$

$$+p_{k}\sigma^{2}_{k}\sum_{m=1}^{m_{k}}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}\bigg{(}\frac{u_{k,i}(1)}{\mu_{k}}+\frac{u_{k,i}(0)}{1-\mu_{k}}\bigg{)}\Bigg{)}^{2}.$$

We will next show that

$$z^{\prime}H_{k}\frac{\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\Sigma$\kern-0.18004pt}}}_{k}}{n_{k}p_{k}q_{k}v_{k}}H_{k}z-\frac{v_{k}^{\rm cluster}}{v_{k}}\stackrel{{\scriptstyle p}}{{\longrightarrow}}0.$$

Given the $\mu_{k}(1-\mu_{k})$ is bounded away from zero, by the weak law of large numbers for arrays, it is enough to show

$$\frac{1}{(n_{k}p_{k}q_{k}v_{k})^{2}}\sum_{m=1}^{m_{k}}E\bigg{[}\bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}(W_{k,i}-\mu_{k})U_{k,i}\bigg{)}^{4}\bigg{]}\longrightarrow 0.$$

Applying the multinomial theorem and the fact that all moments of $W_{k,i}$ as well as all potential outcomes are bounded, we obtain:

$$\frac{1}{(n_{k}p_{k}q_{k}v_{k})^{2}} \sum_{m=1}^{m_{k}}E\bigg{[}\bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}(W_{k,i}-\mu_{k})U_{k,i}\bigg{)}^{4}\bigg{]}$$

$$\leq\frac{c}{(n_{k}p_{k}q_{k}v_{k})^{2}}\Big{(}n_{k}p_{k}q_{k}+n_{k}p_{k}^{2}q_{k}\max_{m}n_{k,m}+n_{k}p_{k}^{3}q_{k}\max_{m}n^{2}_{k,m}+n_{k}p_{k}^{4}q_{k}\max_{m}n^{3}_{k,m}\Big{)}.$$

Now, using $\limsup_{k\rightarrow\infty}\max_{m}n_{k,m}/\min_{m}n_{k,m}<\infty$, $\limsup_{k\rightarrow\infty}p_{k}\min_{m}n_{k,m}/v_{k}<\infty$, and $q_{k}m_{k}\rightarrow\infty$ we obtain

$$\frac{1}{(n_{k}p_{k}q_{k}v_{k})^{2}} \sum_{m=1}^{m_{k}}E\bigg{[}\bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}(W_{k,i}-\mu_{k})U_{k,i}\bigg{)}^{4}\bigg{]}$$

$$\leq c\left(\frac{1}{n_{k}p_{k}q_{k}v_{k}^{2}}+\frac{\max_{m}n_{k,m}}{\min_{m}n_{k,m}}\frac{1}{q_{k}m_{k}v_{k}^{2}}+\frac{p_{k}\max_{m}n^{2}_{k,m}}{v_{k}\min_{m}n_{k,m}}\frac{1}{q_{k}m_{k}v_{k}}+\frac{p_{k}^{2}\max_{m}n^{3}_{k,m}}{v^{2}_{k}\min_{m}n_{k,m}}\frac{1}{q_{k}m_{k}}\right)$$

$$\longrightarrow 0.$$

As a result,

$$\frac{\widehat{V}_{k}^{\rm{cluster}}}{v_{k}}=\frac{v_{k}^{\rm cluster}}{v_{k}}+\scaleto{\mathcal{O}}{5pt}_{p}(1).$$

The robust (sandwich) estimator of the variance of $\sqrt{N_{k}}(\widehat{\tau}_{k}-\tau_{k})$ is given by

$$\widehat{V}_{k}^{\rm robust}=N_{k}z^{\prime}\widehat{Q}_{k}^{-1}\widehat{\Omega}_{k}\widehat{Q}_{k}^{-1}z.$$

where

$$\widehat{\Omega}_{k}=\sum_{i=1}^{n_{k}}R_{k,i}\left(\begin{array}[]{c}\widehat{U}_{k,i}\\ W_{k,i}\widehat{U}_{k,i}\end{array}\right)\left(\begin{array}[]{c}\widehat{U}_{k,i}\\ W_{k,i}\widehat{U}_{k,i}\end{array}\right)^{\prime}.$$

We will derive the limit of $\widehat{V}_{k}^{\rm robust}/v_{k}$. Let

$$\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\Omega$\kern-0.18004pt}}}_{k}=\sum_{i=1}^{n_{k}}R_{k,i}\left(\begin{array}[]{c}U_{k,i}\\ W_{k,i}U_{k,i}\end{array}\right)\left(\begin{array}[]{c}U_{k,i}\\ W_{k,i}U_{k,i}\end{array}\right)^{\prime}.$$

Because potential outcomes (and $W_{k,i}$) are bounded, we obtain

$$\frac{1}{n_{k}p_{k}q_{k}v_{k}}\|\widehat{\Omega}_{k}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\Omega$\kern-0.18004pt}}}_{k}\|\leq c\left(\frac{1}{n_{k}p_{k}q_{k}v_{k}}\sum_{i=1}^{n_{k}}R_{k,i}\right)\max_{i=1,\ldots,n_{k}}|\widehat{U}_{k,i}^{2}-U_{k,i}^{2}|.$$

Because the limsup of the expectation of the first factor (which is non-negative) is bounded and the second factor converges to zero in probability as proved above, we obtain

$$\frac{1}{n_{k}p_{k}q_{k}v_{k}}\|\widehat{\Omega}_{k}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\Omega$\kern-0.18004pt}}}_{k}\|=\scaleto{\mathcal{O}}{5pt}_{p}(1).$$

Notice that

$$\frac{1}{n_{k}p_{k}q_{k}v_{k}}\|\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\Omega$\kern-0.18004pt}}}_{k}\| \leq c\left(\frac{1}{n_{k}p_{k}q_{k}v_{k}}\sum_{i=1}^{n_{k}}R_{k,i}\right)$$

Again, the limsup of the expectation of the right-hand side of this equation is non-negative and bounded. As a result, we obtain $\|\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\Omega$\kern-0.18004pt}}}_{k}\|/(n_{k}p_{k}q_{k}v_{k})=\mathcal{O}_{p}(1)$.

$$\widehat{V}_{k}^{\rm{robust}}/v_{k} =z^{\prime}H_{k}\frac{\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\Omega$\kern-0.18004pt}}}_{k}}{n_{k}p_{k}q_{k}v_{k}}H_{k}z+\scaleto{\mathcal{O}}{5pt}_{p}(1)$$

$$=\frac{1}{n_{k}p_{k}q_{k}v_{k}}\left(\frac{1}{\mu_{k}(1-\mu_{k})}\right)^{2}\sum_{i=1}^{n_{k}}R_{k,i}(W_{k,i}-\mu_{k})^{2}U_{k,i}^{2}+\scaleto{\mathcal{O}}{5pt}_{p}(1).$$

Notice that

$$E\bigg{[}\sum_{i=1}^{n_{k}} R_{k,i}(W_{k,i}-\mu_{k})^{2}U_{k,i}^{2}\bigg{]}=\sum_{i=1}^{n_{k}}p_{k}q_{k}\mu_{k}(1-\mu_{k})\Big{(}(1-\mu_{k})u^{2}_{k,i}(1)+\mu_{k}u^{2}_{k,i}(0)\Big{)}.$$

Finally, notice that

$$\frac{1}{(n_{k}p_{k}q_{k}v_{k})^{2}}\sum_{m=1}^{m_{k}}E\Bigg{[}\Bigg{(}\sum_{i=1}^{n_{k}}R_{k,i}(W_{k,i}-\mu_{k})^{2}U_{k,i}^{2}\Bigg{)}^{2}\Bigg{]} \leq c\frac{n_{k}p_{k}q_{k}+n_{k}p_{k}^{2}q_{k}\max_{m}n_{k,m}}{(n_{k}p_{k}q_{k}v_{k})^{2}}$$

$$\leq c\Bigg{(}\frac{1}{n_{k}p_{k}q_{k}v_{k}^{2}}+\frac{\max_{m}n_{k,m}}{\min_{m}n_{k,m}}\frac{1}{q_{k}m_{k}v_{k}^{2}}\Bigg{)}$$

$$\longrightarrow 0.$$

Therefore, by the weak law of large numbers for arrays, we obtain

$$\frac{\widehat{V}_{k}^{\rm{robust}}}{v_{k}}=\frac{v_{k}^{\rm robust}}{v_{k}}+\scaleto{\mathcal{O}}{5pt}_{p}(1),$$

where

$$v_{k}^{\rm robust}=\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}\bigg{(}\frac{u^{2}_{k,i}(1)}{\mu_{k}}+\frac{u^{2}_{k,i}(0)}{1-\mu_{k}}\bigg{)}.$$

A.3 Fixed effects

A.3.1 Large $k$ distribution

Let

$$\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$N$\kern-0.18004pt}}}_{k,m}=\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}$$

and

$$\widehat{\tau}_{k}^{\rm{\,fixed}} =\frac{\displaystyle\sum_{m=1}^{m_{k}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}Y_{k,i}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})}{\displaystyle\sum_{m=1}^{m_{k}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})},$$ (A.5)

where

$$\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}=\frac{1}{\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$N$\kern-0.18004pt}}}_{k,m}\vee 1}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}.$$

Notice that we need $\liminf_{k\rightarrow\infty}\mu_{k}(1-\mu_{k})-\sigma^{2}_{k}=\liminf_{k\rightarrow\infty}E[A_{k,m}(1-A_{k,m})]>0$ for this estimator to be well-defined in large samples (otherwise, the denominator in the formula for $\widehat{\tau}_{k}^{\rm{\,fixed}}$ could be equal to zero). Although it is not strictly necessary, and because it entails little loss of generality and simplifies the exposition, we will assume that the supports of the cluster probabilities, $A_{k,m}$, are bounded away from zero and one (uniformly in $k$ and $m$). In finite samples we assign $\widehat{\tau}_{k}^{\rm{\,fixed}}=0$ to the cases when the denominator of $\widehat{\tau}_{k}^{\rm{\,fixed}}$ in equation (A.5) is equal to zero. Notice that

$$\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})=0.$$

Let

$$\alpha_{k,m}=\frac{1}{n_{k,m}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}y_{k,i}(0),\quad\tau_{k,m}=\frac{1}{n_{k,m}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}(y_{k,i}(1)-y_{k,i}(0)),$$

$e_{k,i}(0)=y_{k,i}(0)-\alpha_{k,m_{k,i}}$, and $e_{k,i}(1)=y_{k,i}(1)-\alpha_{k,m_{k,i}}-\tau_{k,m_{k,i}}$. It follows that

$$\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}e_{k,i}(1)=\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}e_{k,i}(0)=0.$$

Now, $Y_{k,i}=e_{k,i}(1)W_{k,i}+e_{k,i}(0)(1-W_{k,i})+\alpha_{k,m_{k,i}}+\tau_{k,m_{k,i}}W_{k,i}$. Then,

$$\widehat{\tau}_{k}^{\rm{\,fixed}} =\frac{\displaystyle\sum_{m=1}^{m_{k}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}((e_{k,i}(1)+\tau_{k,m})W_{k,i}+e_{k,i}(0)(1-W_{k,i}))(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})}{\displaystyle\sum_{m=1}^{m_{k}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})}.$$

Let

$$\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\tau$\kern-0.18004pt}}}_{k} =\frac{\displaystyle\sum_{m=1}^{m_{k}}\tau_{k,m}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})}{\displaystyle\sum_{m=1}^{m_{k}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})},$$ (A.6)

where, as before, we make $\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\tau$\kern-0.18004pt}}}_{k}=0$ if the denominator on the right-hand side of (A.6) is equal to zero. Now, $\widehat{\tau}_{k}^{\rm{\,fixed}}-\tau_{k}=(\widehat{\tau}_{k}^{\rm{\,fixed}}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\tau$\kern-0.18004pt}}}_{k})+(\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\tau$\kern-0.18004pt}}}_{k}-\tau_{k})$, where

$$\widehat{\tau}_{k}^{\rm{\,fixed}}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\tau$\kern-0.18004pt}}}_{k}=\frac{\displaystyle\sum_{m=1}^{m_{k}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}(e_{k,i}(1)W_{k,i}+e_{k,i}(0)(1-W_{k,i}))(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})}{\displaystyle\sum_{m=1}^{m_{k}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})}$$

and

$$\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\tau$\kern-0.18004pt}}}_{k}-\tau_{k}=\frac{\displaystyle\sum_{m=1}^{m_{k}}(\tau_{k,m}-\tau_{k})\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})}{\displaystyle\sum_{m=1}^{m_{k}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})}.$$

Notice that outcomes enter the term $\widehat{\tau}_{k}^{\rm{\,fixed}}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\tau$\kern-0.18004pt}}}_{k}$ only through the intra-cluster errors, $e_{k,i}(1)$ and $e_{k,i}(0)$. In contrast, the term $\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\tau$\kern-0.18004pt}}}_{k}-\tau_{k}$ depends on outcomes only through inter-cluster variability in treatment effects, $\tau_{k,m}-\tau_{k}$. The numerator in the expression for $\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\tau$\kern-0.18004pt}}}_{k}-\tau_{k}$ in the last displayed equation does not have mean zero in general, and this will be reflected in a bias term, $B_{k}$, which we define next. Let,

$$D_{k}=\frac{1}{n_{k}p_{k}q_{k}}\sum_{m=1}^{m_{k}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}),$$

and

$$B_{k}=-\frac{\displaystyle\frac{1}{n_{k}p_{k}}E[A_{k,m}(1-A_{k,m})]\sum_{m=1}^{m_{k}}(\tau_{k,m}-\tau_{k})(1-(1-p_{k})^{n_{k,m}})}{\displaystyle\frac{1}{n_{k}p_{k}q_{k}}\sum_{m=1}^{m_{k}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})}.$$

Then, $\sqrt{n_{k}p_{k}q_{k}}(\widehat{\tau}_{k}^{\rm{\,fixed}}-\tau_{k}-B_{k})=F_{k}/D_{k}$, where

$$F_{k}=\sum_{m=1}^{m_{k}}(\psi_{k,m}-\hbox{\vbox{\hrule height=0.5pt\kern 0.77498pt\hbox{\kern-0.18004pt$\psi$\kern-0.18004pt}}}_{k,m})+(\varphi_{k,m}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\varphi$\kern-0.18004pt}}}_{k,m}),$$

$$\psi_{k,m} =\frac{1}{\sqrt{n_{k}p_{k}q_{k}}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}(e_{k,i}(1)W_{k,i}+e_{k,i}(0)(1-W_{k,i}))(W_{k,i}-A_{k,m}),$$

$$\hbox{\vbox{\hrule height=0.5pt\kern 0.77498pt\hbox{\kern-0.18004pt$\psi$\kern-0.18004pt}}}_{k,m} =\frac{1}{\sqrt{n_{k}p_{k}q_{k}}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}(e_{k,i}(1)W_{k,i}+e_{k,i}(0)(1-W_{k,i}))(\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}-A_{k,m}),$$

$$\varphi_{k,m} =\frac{1}{\sqrt{n_{k}p_{k}q_{k}}}(\tau_{k,m}-\tau_{k})\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}\big{(}R_{k,i}W_{k,i}(W_{k,i}-A_{k,m})-p_{k}q_{k}E[A_{k,m}(1-A_{k,m})]\big{)},$$

and

$$\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\varphi$\kern-0.18004pt}}}_{k,m} =\frac{1}{\sqrt{n_{k}p_{k}q_{k}}}(\tau_{k,m}-\tau_{k})\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}(\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}-A_{k,m})$$

$$\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad-q_{k}E[A_{k,m}(1-A_{k,m})](1-(1-p_{k})^{n_{k,m}})\Bigg{)}.$$

The terms $\psi_{k,m}$ and $\hbox{\vbox{\hrule height=0.5pt\kern 0.77498pt\hbox{\kern-0.18004pt$\psi$\kern-0.18004pt}}}_{k,m}$ depend on the within-cluster errors $e_{k,i}(1)$ and $e_{k,i}(0)$. The terms $\varphi_{k,m}$ and $\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\varphi$\kern-0.18004pt}}}_{k,m}$ depend on the inter-clusters errors $\tau_{k,m}-\tau_{k}$. $\psi_{k,m}$ and $\varphi_{k,m}$ replace $\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}$ with $A_{k,m}$, while $\hbox{\vbox{\hrule height=0.5pt\kern 0.77498pt\hbox{\kern-0.18004pt$\psi$\kern-0.18004pt}}}_{k,m}$ and $\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\varphi$\kern-0.18004pt}}}_{k,m}$ correct for the difference, $\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}-A_{k,m}$.

It can be seen (in intermediate calculations below) that

$$E\Bigg{[}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}(W_{k,i}-A_{k,m})\Bigg{]}=n_{k,m}p_{k}q_{k}E[A_{k,m}(1-A_{k,m})]$$

and

$$E\Bigg{[}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}(\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}-A_{k,m})\Bigg{]} =q_{k}E[A_{k,m}(1-A_{k,m})](1-(1-p_{k})^{n_{k,m}}).$$

These two expectations are substracted in $\varphi_{k,m}$ and $\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\varphi$\kern-0.18004pt}}}_{k,m}$ , so $\varphi_{k,m}$ and $\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\varphi$\kern-0.18004pt}}}_{k,m}$ have mean zero. Doing so for $\varphi_{k,m}$ does not require adjustments elsewhere. Because

$$\sum_{m=1}^{m_{k}}(\tau_{k,m}-\tau_{k})n_{k,m}=0,$$

the $n_{k,m}p_{k}q_{k}E[A_{k,m}(1-A_{k,m})]$ terms do not change the sum $F_{k}$. In contrast, demeaning $\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\varphi$\kern-0.18004pt}}}_{k,m}$ creates the bias term $B_{k}$. If the size of the clusters $n_{k,m}$ does not vary across clusters, then $B_{k}$ is equal to zero. More generally, $\sqrt{n_{k}p_{k}q_{k}}D_{k}B_{k}=\mathcal{O}(m_{k}\sqrt{q_{k}/(n_{k}p_{k})})$. Therefore, if

$$\frac{m_{k}q_{k}}{p_{k}(n_{k}/m_{k})}\longrightarrow 0,$$ (A.7)

(that is, if the expected number of sampled clusters is small relative to the expected number of sampled observations per sampled cluster) then $\sqrt{n_{k}p_{k}q_{k}}D_{k}B_{k}$ converges to zero. As a result, $\sqrt{n_{k}p_{k}q_{k}}B_{k}$ converges in probability to zero, because, as we will show later, $D_{k}$ converges in probability to $\mu_{k}(1-\mu_{k})-\sigma_{k}^{2}$, which is bounded away from zero. In our large sample analysis, we will assume that the expected number of sampled clusters grows to infinity, $m_{k}q_{k}\rightarrow\infty$. Then, equation (A.7) implies that the expected number of observations per sampled cluster goes to infinity, $p_{k}(n_{k}/m_{k})\rightarrow\infty$. Notice also that $n_{k}p_{k}q_{k}=(n_{k}p_{k}/m_{k})(m_{k}q_{k})\rightarrow\infty$.

We summarize now the assumptions we made thus far. We first assumed that the supports of the cluster probabilities, $A_{k,m}$, are bounded away from zero and one (uniformly in $k$ and $m$), and that potential outcomes are bounded. Moreover, we assumed $m_{k}q_{k}\rightarrow\infty$ and $(m_{k}q_{k})/((p_{k}n_{k})/m_{k})\rightarrow 0$. These imply $(p_{k}n_{k})/m_{k}\rightarrow\infty$ and $n_{k}p_{k}q_{k}\rightarrow\infty$. We will add the assumption that the ratio between maximum and minimum cluster size is bounded, $\limsup_{k\rightarrow\infty}\max_{m}n_{k,m}/\min_{m}n_{k,m}<\infty$. This assumption implies $p_{k}\min_{m}n_{k,m}\rightarrow\infty$ and $(m_{k}q_{k})/(p_{k}\min_{m}n_{k,m})\rightarrow 0$.

We will now study the behavior of $D_{k}$. Notice that

$$E\Bigg{[}\sum_{m=1}^{m_{k}}\sum_{i=1}^{n_{k}} 1\{m_{k,i}=m\}R_{k,i}W_{k,i}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})\Bigg{]}$$

$$=E\Bigg{[}\sum_{m=1}^{m_{k}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}(W_{k,i}-A_{k,m})\Bigg{]}$$

$$-E\Bigg{[}\sum_{m=1}^{m_{k}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}(\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}-A_{k,m})\Bigg{]}$$

$$=n_{k}p_{k}q_{k}E[A_{k,m}(1-A_{k,m})]-q_{k}E[A_{k,m}(1-A_{k,m})]\sum_{m=1}^{m_{k}}(1-(1-p_{k})^{n_{k,m}}).$$

In addition,

$$\frac{1}{(n_{k}p_{k}q_{k})^{2}}\sum_{m=1}^{m_{k}}E\Bigg{[}\Bigg{(} \sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})\Bigg{)}^{2}\Bigg{]}$$

$$\leq c\,\frac{n_{k}p_{k}q_{k}+n_{k}p_{k}^{2}q_{k}\max_{m}n_{k,m}}{(n_{k}p_{k}q_{k})^{2}}$$

$$=c\,\Bigg{(}\frac{1}{n_{k}p_{k}q_{k}}+\frac{\max_{m}n_{k,m}}{\min_{m}n_{k,m}}\frac{1}{m_{k}q_{k}}\Bigg{)}\longrightarrow 0.$$

The weak law of large numbers for arrays implies

$$D_{k}-E[A_{k,m}(1-A_{k,m})]+\frac{1}{n_{k}p_{k}}E[A_{k,m}(1-A_{k,m})]\sum_{m=1}^{m_{k}}(1-(1-p_{k})^{n_{k,m}})\stackrel{{\scriptstyle p}}{{\longrightarrow}}0.$$

Because $m_{k}/(n_{k}p_{k})\rightarrow 0$ and $E[A_{k,m}(1-A_{k,m})]=\mu_{k}(1-\mu_{k})-\sigma_{k}^{2}$, we obtain

$$D_{k}-(\mu_{k}(1-\mu_{k})-\sigma_{k}^{2})\stackrel{{\scriptstyle p}}{{\longrightarrow}}0.$$

We now turn our attention to $F_{k}$. We will first calculate the variance of $\psi_{k,m}$. Let $Q_{k,m}$ be a binary variable that takes value one if cluster $m$ in population $k$ is sampled, and zero otherwise. Notice that

$$E[R_{k,i}W_{k,i}(W_{k,i}-A_{k,m})|A_{k,m},Q_{k,m}=1,m_{k,i}=m]=p_{k}A_{k,m}(1-A_{k,m}),$$

and

$$E[R_{k,i}(1-W_{k,i})(W_{k,i}-A_{k,m})|A_{k,m},Q_{k,m}=1,m_{k,i}=m]=-p_{k}A_{k,m}(1-A_{k,m}).$$

Consider now

$$\psi_{k,m,1} =\frac{1}{\sqrt{n_{k}p_{k}q_{k}}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}(W_{k,i}-A_{k,m})e_{k,i}(1)$$

$$=\frac{Q_{k,m}}{\sqrt{n_{k}p_{k}q_{k}}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}\Big{(}R_{k,i}W_{k,i}(W_{k,i}-A_{k,m})-p_{k}A_{k,m}(1-A_{k,m})\Big{)}e_{k,i}(1),$$

and

$$\psi_{k,m,0} =\frac{1}{\sqrt{n_{k}p_{k}q_{k}}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}(1-W_{k,i})(W_{k,i}-A_{k,m})e_{k,i}(0)$$

$$=\frac{Q_{k,m}}{\sqrt{n_{k}p_{k}q_{k}}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}\Big{(}R_{k,i}(1-W_{k,i})(W_{k,i}-A_{k,m})+p_{k}A_{k,m}(1-A_{k,m})\Big{)}e_{k,i}(0).$$

It holds that $\psi_{k,m}=\psi_{k,m,1}+\psi_{k,m,0}$ and $E[\psi_{k,m}]=0$. Now, notice that

$$E[\psi_{k,m,1}^{2}] =\frac{1}{n_{k}}E[A_{k,m}(1-A_{k,m})^{2}-p_{k}A_{k,m}^{2}(1-A_{k,m})^{2}]\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}e_{k,i}^{2}(1),$$

$$E[\psi_{k,m,0}^{2}] =\frac{1}{n_{k}}E[A^{2}_{k,m}(1-A_{k,m})-p_{k}A_{k,m}^{2}(1-A_{k,m})^{2}]\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}e_{k,i}^{2}(0),$$

and

$$E[\psi_{k,m,1}\psi_{k,m,0}] =\frac{1}{n_{k}}p_{k}E[A_{k,m}^{2}(1-A_{k,m})^{2}]\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}e_{k,i}(1)e_{k,i}(0).$$

Therefore,

$$E[(\psi_{k,m,1}+\psi_{k,m,0})^{2}] =\frac{1}{n_{k}}E[A_{k,m}(1-A_{k,m})^{2}]\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}e_{k,i}^{2}(1)$$

$$+\frac{1}{n_{k}}E[A^{2}_{k,m}(1-A_{k,m})]\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}e_{k,i}^{2}(0)$$

$$-\frac{1}{n_{k}}p_{k}E[A^{2}_{k,m}(1-A_{k,m})^{2}]\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}(e_{k,i}(1)-e_{k,i}(0))^{2},$$

and

$$\sum_{m=1}^{m_{k}}E[(\psi_{k,m,1}+\psi_{k,m,0})^{2}] =E[A_{k,m}(1-A_{k,m})^{2}]\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}e_{k,i}^{2}(1)$$

$$+E[A^{2}_{k,m}(1-A_{k,m})]\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}e_{k,i}^{2}(0)$$

$$-p_{k}E[A^{2}_{k,m}(1-A_{k,m})^{2}]\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}(e_{k,i}(1)-e_{k,i}(0))^{2}.$$ (A.8)

We will next show that the terms $\hbox{\vbox{\hrule height=0.5pt\kern 0.77498pt\hbox{\kern-0.18004pt$\psi$\kern-0.18004pt}}}_{k,m}$ do not matter for the asymptotic distribution of $\sqrt{n_{k}p_{k}q_{k}}(\widehat{\tau}_{k}-\tau_{k})$. Notice that, because the cluster sum of $e_{k,i}(1)$ is equal to zero, we obtain $E[\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\psi$\kern-0.18004pt}}}_{k,m}]=0$ and, therefore,

$$\sum_{m=1}^{m_{k}}E\Big{[}\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\psi$\kern-0.18004pt}}}_{k,m}\Big{]}=0.$$

Moreover

$$2\sum_{i=1}^{n_{k}-1}\sum_{j=i+1}^{n_{k}}1\{m_{k,i}=m_{k,j}=m\}e_{k,i}(1)e_{k,j}(1)=-\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}e_{k,i}^{2}(1)\leq 0.$$

In addition, $E[R_{k,i}W_{k,i}(\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}-A_{k,m})^{2}|m_{k,i}=m]\leq q_{k}E[A_{k,m}(1-A_{k,m})]/n_{k,m}$ (see intermediate calculations). Therefore,

$$E\Bigg{[}\Bigg{(} \sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}(\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}-A_{k,m})e_{k,i}(1)\Bigg{)}^{2}\Bigg{]}$$

$$=\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}E\Big{[}R_{k,i}W_{k,i}(\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}-A_{k,m})^{2}|m_{k,i}=m\Big{]}e_{k,i}^{2}(1)$$

$$+2\sum_{i=1}^{n_{k}-1}\sum_{j=i+1}^{n_{k}}E\Big{[}1\{m_{k,i}=m_{k,j}=m\}R_{k,i}R_{k,j}W_{k,i}W_{k,j}(\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}-A_{k,m})^{2}\Big{]}e_{k,i}(1)e_{k,j}(1)$$

$$\leq q_{k}E[A_{k,m}(1-A_{k,m})]\frac{1}{n_{k,m}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}e_{k,i}^{2}(1).$$

Now, because errors are bounded, we obtain

$$\sum_{m=1}^{m_{k}}E\Bigg{[}\Bigg{(}\frac{1}{\sqrt{n_{k}p_{k}q_{k}}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}(\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}-A_{k,m})e_{k,i}(1)\Bigg{)}^{2}\Bigg{]}\leq c\,\frac{m_{k}}{n_{k}p_{k}}.$$ (A.9)

Because $m_{k}/(n_{k}p_{k})\rightarrow 0$, the weak law of large numbers for arrays, implies,

$$\frac{1}{\sqrt{n_{k}p_{k}q_{k}}}\sum_{m=1}^{m_{k}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}(\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}-A_{k,m})e_{k,i}(1)\stackrel{{\scriptstyle p}}{{\longrightarrow}}0.$$

with the analogous result involving the errors $e_{k,i}(0)$. If follows that

$$\sum_{m=1}^{m_{k}}\hbox{\vbox{\hrule height=0.5pt\kern 0.77498pt\hbox{\kern-0.18004pt$\psi$\kern-0.18004pt}}}_{k,m}\stackrel{{\scriptstyle p}}{{\longrightarrow}}0.$$

Consider now $\varphi_{k,m}$. Notice that

$$E\Big{[}\Big{(}R_{k,i}W_{k,i} (W_{k,i}-A_{k,m})-p_{k}q_{k}E[A_{k,m}(1-A_{k,m})])\Big{)}^{2}\Big{]}$$

$$=p_{k}q_{k}E[A_{k,m}(1-A_{k,m})^{2}]-p_{k}^{2}q_{k}^{2}\big{(}E[A_{k,m}(1-A_{k,m})]\big{)}^{2},$$

and

$$E\Big{[}\Big{(}R_{k,i}W_{k,i} (W_{k,i}-A_{k,m})-p_{k}q_{k}E[A_{k,m}(1-A_{k,m})])\Big{)}$$

$$\times\Big{(}R_{k,j}W_{k,j}(W_{k,j}-A_{k,m})-p_{k}q_{k}E[A_{k,m}(1-A_{k,m})])\Big{)}\big{|}m_{k,i}=m_{k,j}=m\Big{]}$$

$$\qquad=p_{k}^{2}q_{k}E[A^{2}_{k,m}(1-A_{k,m})^{2}]-p_{k}^{2}q_{k}^{2}\big{(}E[A_{k,m}(1-A_{k,m})]\big{)}^{2}.$$

Therefore,

$$E[\varphi_{k,m}^{2}] =\Big{(}E[A_{k,m}(1-A_{k,m})^{2}]-p_{k}q_{k}(E[A_{k,m}(1-A_{k,m})])^{2}\Big{)}\frac{n_{k,m}}{n_{k}}(\tau_{k,m}-\tau_{k})^{2}$$

$$+\Big{(}p_{k}E[A^{2}_{k,m}(1-A_{k,m})^{2}]-p_{k}q_{k}(E[A_{k,m}(1-A_{k,m})])^{2}\Big{)}\frac{n_{k,m}(n_{k,m}-1)}{n_{k}}(\tau_{k,m}-\tau_{k})^{2},$$

and

$$\sum_{m=1}^{m_{k}}E[\varphi_{k,m}^{2}] =\Big{(}E[A_{k,m}(1-A_{k,m})^{2}]-p_{k}q_{k}(E[A_{k,m}(1-A_{k,m})])^{2}\Big{)}\sum_{m=1}^{m_{k}}\frac{n_{k,m}}{n_{k}}(\tau_{k,m}-\tau_{k})^{2}$$

$$+\Big{(}p_{k}E[A^{2}_{k,m}(1-A_{k,m})^{2}]-p_{k}q_{k}(E[A_{k,m}(1-A_{k,m})])^{2}\Big{)}\sum_{m=1}^{m_{k}}\frac{n_{k,m}(n_{k,m}-1)}{n_{k}}(\tau_{k,m}-\tau_{k})^{2}.$$

Next, we calculate the variance of $\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\varphi$\kern-0.18004pt}}}_{k,m}$. Using results on the moments of a Binomial distribution, we obtain, for $n\geq 1$,

$$E\Bigg{[}\Bigg{(}\sum_{i=1}^{n_{k}} 1\{m_{k,i}=m\}R_{k,i}W_{k,i}(\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}-A_{k,m})\Bigg{)}^{2}\Big{|}\begin{array}[]{l}Q_{k,m}=1,\\ \hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$N$\kern-0.18004pt}}}_{k,m}=n\end{array}\Bigg{]}$$

$$=\frac{1}{n^{2}}E\Bigg{[}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}\Big{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}-nA_{k,m}\Big{)}\Bigg{)}^{2}\Big{|}\begin{array}[]{l}Q_{k,m}=1,\\ \hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$N$\kern-0.18004pt}}}_{k,m}=n\end{array}\Bigg{]}$$

$$=nE[A_{k,m}^{3}(1-A_{k,m})]+E[A^{2}_{k,m}(1-A_{k,m})(5-7A_{k,m})]$$

$$+\frac{1}{n}E[A_{k,m}(1-A_{k,m})(6A^{2}_{k,m}-6A_{k,m}+1)].$$

Therefore,

$$E\Bigg{[}\Bigg{(}\sum_{i=1}^{n_{k}} 1\{m_{k,i}=m\}R_{k,i}W_{k,i}(\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}-A_{k,m})\Bigg{)}^{2}\Bigg{]}$$

$$=n_{k,m}p_{k}q_{k}E[A_{k,m}^{3}(1-A_{k,m})]+q_{k}E[A^{2}_{k,m}(1-A_{k,m})(5-7A_{k,m})](1-(1-p_{k})^{n_{k,m}})$$

$$+q_{k}E[A_{k,m}(1-A_{k,m})(6A^{2}_{k,m}-6A_{k,m}+1)]r_{k,m},$$

where

$$r_{k,m}=\sum_{n=1}^{n_{k,m}}\frac{1}{n}\Pr(\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$N$\kern-0.18004pt}}}_{k,m}=n|Q_{k,m}=1)\leq\sum_{n=1}^{n_{k,m}}\Pr(\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$N$\kern-0.18004pt}}}_{k,m}=n|Q_{k,m}=1)\leq 1.$$

It follows that,

$$E[\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\varphi$\kern-0.18004pt}}}^{2}_{k,m}] =(\tau_{k,m}-\tau_{k})^{2}\Big{(}\frac{n_{k,m}}{n_{k}}E[A_{k,m}^{3}(1-A_{k,m})]+\frac{1}{n_{k}p_{k}}E[A^{2}_{k,m}(1-A_{k,m})(5-7A_{k,m})](1-(1-p_{k})^{n_{k,m}})$$

$$+\frac{1}{n_{k}p_{k}}E[A_{k,m}(1-A_{k,m})(6A^{2}_{k,m}-6A_{k,m}+1)]r_{k,m}$$

$$-\frac{q_{k}}{n_{k}p_{k}}(E[A_{k,m}(1-A_{k,m})])^{2}(1-(1-p_{k})^{n_{k,m}})^{2}\Big{)}.$$

Therefore,

$$\sum_{m=1}^{m_{k}}E[\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\varphi$\kern-0.18004pt}}}^{2}_{k,m}]=\sum_{m=1}^{m_{k}}(\tau_{k,m}-\tau_{k})^{2}\Big{(}\frac{n_{k,m}}{n_{k}}\Big{)}E[A_{k,m}^{3}(1-A_{k,m})]+\scaleto{\mathcal{O}}{5pt}(1).$$

We will now study the covariance between $\varphi_{k,m}$ and $\hbox{\vbox{\hrule height=0.5pt\kern 0.77498pt\hbox{\kern-0.18004pt$\varphi$\kern-0.18004pt}}}_{k,m}$. Using results on the moments of a Binomial distribution, we obtain, for $n\geq 1$,

$$E\Bigg{[}\Bigg{(} \sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}(W_{k,i}-A_{k,m})\Bigg{)}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}(\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}-A_{k,m})\Bigg{)}\Big{|}\begin{array}[]{l}Q_{k,m}=1,\\ \hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$N$\kern-0.18004pt}}}_{k,m}=n\end{array}\Bigg{]}$$

$$=E\Bigg{[}\frac{1-A_{k,m}}{n}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}\Bigg{)}^{2}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}-nA_{k,m})\Bigg{)}\Big{|}\begin{array}[]{l}Q_{k,m}=1,\\ \hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$N$\kern-0.18004pt}}}_{k,m}=n\end{array}\Bigg{]}$$

$$=2nE[A^{2}_{k,m}(1-A_{k,m})^{2}]+E[A_{k,m}(1-A_{k,m})^{2}(1-2A_{k,m})].$$

Therefore,

$$E\Bigg{[}\Bigg{(} \sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}(W_{k,i}-A_{k,m})\Bigg{)}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}(\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}-A_{k,m})\Bigg{)}\Bigg{]}$$

$$=2n_{k,m}p_{k}q_{k}E[A^{2}_{k,m}(1-A_{k,m})^{2}]+q_{k}E[A_{k,m}(1-A_{k,m})^{2}(1-2A_{k,m})]\Pr(\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$N$\kern-0.18004pt}}}_{k,m}\geq 1|Q_{k,m}=1).$$

In addition,

$$E\Bigg{[} \sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}(W_{k,i}-A_{k,m})\Bigg{]}E\Bigg{[}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}(\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}-A_{k,m})\Bigg{]}$$

$$=n_{k,m}p_{k}q^{2}_{k}(E[A_{k,m}(1-A_{k,m})])^{2}\Pr(\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$N$\kern-0.18004pt}}}_{k,m}\geq 1|Q_{k,m}=1).$$

As a result,

$$E[\varphi_{k,m}\hbox{\vbox{\hrule height=0.5pt\kern 0.77498pt\hbox{\kern-0.18004pt$\varphi$\kern-0.18004pt}}}_{k,m}] =\Big{(}2E[A^{2}_{k,m}(1-A_{k,m})^{2}]-q_{k}(E[A_{k,m}(1-A_{k,m})])^{2}\Big{)}(\tau_{k,m}-\tau_{k})^{2}\left(\frac{n_{k,m}}{n_{k}}\right)+\mathcal{O}\left(\frac{1}{n_{k}p_{k}}\right)$$

$$+\mathcal{O}\left(\frac{q_{k}}{p_{k}\min_{m}n_{k,m}}\big{(}p_{k}\min_{m}n_{k,m}(1-p_{k})^{\min_{m}n_{k,m}}\big{)}\right).$$

Notice that $m_{k}/(n_{k}p_{k})\rightarrow 0$. In addition, $m_{k}q_{k}/(p_{k}\min_{m}n_{k,m})\rightarrow 0$ and

$$p_{k}\min_{m}n_{k,m}(1-p_{k})^{\min_{m}n_{k,m}} =p_{k}\min_{m}n_{k,m}\left(1-\frac{p_{k}\min_{m}n_{k,m}}{\min_{m}n_{k,m}}\right)^{\min_{m}n_{k,m}}$$

$$<p_{k}\min_{m}n_{k,m}e^{-p_{k}\min_{m}n_{k,m}}\longrightarrow 0.$$

Therefore,

$$\sum_{m=1}^{m_{k}}E[\varphi_{k,m}\hbox{\vbox{\hrule height=0.5pt\kern 0.77498pt\hbox{\kern-0.18004pt$\varphi$\kern-0.18004pt}}}_{k,m}]=\Big{(}2E[A^{2}_{k,m}(1-A_{k,m})^{2}]-q_{k}(E[A_{k,m}(1-A_{k,m})])^{2}\Big{)}\sum_{m=1}^{m_{k}}(\tau_{k,m}-\tau_{k})^{2}\left(\frac{n_{k,m}}{n_{k}}\right)+\ \scaleto{\mathcal{O}}{5pt}(1).$$

Next, we will study the remaining covariances between $\psi_{k,m}$, $\varphi_{k,m}$, $\hbox{\vbox{\hrule height=0.5pt\kern 0.77498pt\hbox{\kern-0.18004pt$\psi$\kern-0.18004pt}}}_{k,m}$, and $\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\varphi$\kern-0.18004pt}}}_{k,m}$. Because the intra-cluster errors, $e_{k,i}(1)$ and $e_{k,i}(0)$ sum to zero, it can be easily seen that $E[\psi_{k,m}\varphi_{k,m}]=E[\psi_{k,m}\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\varphi$\kern-0.18004pt}}}_{k,m}]=0$. It can also be seen that the inter-clusters sums of covariances between $\hbox{\vbox{\hrule height=0.5pt\kern 0.77498pt\hbox{\kern-0.18004pt$\psi$\kern-0.18004pt}}}_{k,m}$ and any of the other terms go to zero. To prove this for the covariance with $\psi_{k,m}$, we have

$$\left(\sum_{m=1}^{m_{k}}E[|\psi_{k,m}\hbox{\vbox{\hrule height=0.5pt\kern 0.77498pt\hbox{\kern-0.18004pt$\psi$\kern-0.18004pt}}}_{k,m}|]\right)^{2} \leq\left(\sum_{m=1}^{m_{k}}(E[\psi^{2}_{k,m}]E[\hbox{\vbox{\hrule height=0.5pt\kern 0.77498pt\hbox{\kern-0.18004pt$\psi$\kern-0.18004pt}}}^{2}_{k,m}])^{1/2}\right)^{2}$$

$$\leq\sum_{m=1}^{m_{k}}E[\psi^{2}_{k,m}]\sum_{m=1}^{m_{k}}E[\hbox{\vbox{\hrule height=0.5pt\kern 0.77498pt\hbox{\kern-0.18004pt$\psi$\kern-0.18004pt}}}^{2}_{k,m}]$$

$$=\mathcal{O}(1)\scaleto{\mathcal{O}}{5pt}(1)=\scaleto{\mathcal{O}}{5pt}(1).$$

The same argument and result applies to $E[\hbox{\vbox{\hrule height=0.5pt\kern 0.77498pt\hbox{\kern-0.18004pt$\psi$\kern-0.18004pt}}}_{k,m}\varphi_{k,m}]$ and $E[\hbox{\vbox{\hrule height=0.5pt\kern 0.77498pt\hbox{\kern-0.18004pt$\psi$\kern-0.18004pt}}}_{k,m}\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\varphi$\kern-0.18004pt}}}_{k,m}]$. Putting all the pieces together, we obtain

$$n_{k}p_{k}q_{k}E[D_{k}^{2}(\widehat{\tau}_{k}^{\rm{\,fixed}}-\tau_{k})^{2}]=f_{k}+\scaleto{\mathcal{O}}{5pt}(1),$$

where

$$f_{k} =E[A_{k,m}(1-A_{k,m})^{2}]\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}e_{k,i}^{2}(1)+E[A^{2}_{k,m}(1-A_{k,m})]\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}e_{k,i}^{2}(0)$$

$$-p_{k}E[A^{2}_{k,m}(1-A_{k,m})^{2}]\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}(e_{k,i}(1)-e_{k,i}(0))^{2}$$

$$+\Big{(}E[A_{k,m}(1-A_{k,m})^{2}]-p_{k}q_{k}(E[A_{k,m}(1-A_{k,m})])^{2}\Big{)}\sum_{m=1}^{m_{k}}\frac{n_{k,m}}{n_{k}}(\tau_{k,m}-\tau_{k})^{2}$$

$$+\Big{(}p_{k}E[A^{2}_{k,m}(1-A_{k,m})^{2}]-p_{k}q_{k}(E[A_{k,m}(1-A_{k,m})])^{2}\Big{)}\sum_{m=1}^{m_{k}}\frac{n_{k,m}(n_{k,m}-1)}{n_{k}}(\tau_{k,m}-\tau_{k})^{2}$$

$$+E[A_{k,m}^{3}(1-A_{k,m})]\sum_{m=1}^{m_{k}}(\tau_{k,m}-\tau_{k})^{2}\Big{(}\frac{n_{k,m}}{n_{k}}\Big{)}$$

$$-2\Big{(}2E[A^{2}_{k,m}(1-A_{k,m})^{2}]-q_{k}(E[A_{k,m}(1-A_{k,m})])^{2}\Big{)}\sum_{m=1}^{m_{k}}(\tau_{k,m}-\tau_{k})^{2}\left(\frac{n_{k,m}}{n_{k}}\right).$$

Collecting terms with identical factors, we obtain

$$f_{k} =E[A_{k,m}(1-A_{k,m})^{2}]\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}e_{k,i}^{2}(1)+E[A^{2}_{k,m}(1-A_{k,m})]\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}e_{k,i}^{2}(0)$$

$$-p_{k}E[A^{2}_{k,m}(1-A_{k,m})^{2}]\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}(e_{k,i}(1)-e_{k,i}(0))^{2}$$

$$+\Big{(}E[A_{k,m}(1-A_{k,m})^{2}]-(4+p_{k})E[A^{2}_{k,m}(1-A_{k,m})^{2}]$$

$$\qquad\qquad+E[A_{k,m}^{3}(1-A_{k,m})]+2q_{k}(E[A_{k,m}(1-A_{k,m})])^{2}\Big{)}\sum_{m=1}^{m_{k}}\frac{n_{k,m}}{n_{k}}(\tau_{k,m}-\tau_{k})^{2}$$

$$+\Big{(}p_{k}E[A^{2}_{k,m}(1-A_{k,m})^{2}]-p_{k}q_{k}(E[A_{k,m}(1-A_{k,m})])^{2}\Big{)}\sum_{m=1}^{m_{k}}\frac{n^{2}_{k,m}}{n_{k}}(\tau_{k,m}-\tau_{k})^{2}.$$

The first three terms in the expression above depend on intra-cluster heterogeneity in potential outcomes and treatment effects. The last two terms depend on inter-cluster variation in average treatment effects.

A more compact expression for $f_{k}$ is

$$f_{k} =E[A_{k,m}(1-A_{k,m})^{2}]\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}e_{k,i}^{2}(1)+E[A^{2}_{k,m}(1-A_{k,m})]\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}e_{k,i}^{2}(0)$$

$$-p_{k}E[A^{2}_{k,m}(1-A_{k,m})^{2}]\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}(e_{k,i}(1)-e_{k,i}(0))^{2}$$

$$+\Big{(}E[A_{k,m}(1-A_{k,m})]-(5+p_{k})E[A^{2}_{k,m}(1-A_{k,m})^{2}]$$

$$\qquad\qquad+2q_{k}(E[A_{k,m}(1-A_{k,m})])^{2}\Big{)}\sum_{m=1}^{m_{k}}\frac{n_{k,m}}{n_{k}}(\tau_{k,m}-\tau_{k})^{2}$$

$$+\Big{(}p_{k}E[A^{2}_{k,m}(1-A_{k,m})^{2}]-p_{k}q_{k}(E[A_{k,m}(1-A_{k,m})])^{2}\Big{)}\sum_{m=1}^{m_{k}}\frac{n^{2}_{k,m}}{n_{k}}(\tau_{k,m}-\tau_{k})^{2}.$$ (A.10)

Notice that the first four terms in (A.3.1) are bounded, and that

$$E[A^{2}_{k,m}(1-A_{k,m})^{2}]-q_{k}(E[A_{k,m}(1-A_{k,m})])^{2}=\mbox{var}(A_{k,m}(1-A_{k,m}))+(1-q_{k})(E[A_{k,m}(1-A_{k,m})])^{2}.$$

Assume that

$$\liminf_{k\rightarrow\infty}\sum_{m=1}^{m_{k}}\frac{n_{k,m}}{n_{k}}(\tau_{k,m}-\tau_{k})^{2}>0,$$ (A.11)

and

$$\liminf_{k\rightarrow\infty}\mbox{var}(A_{k,m}(1-A_{k,m}))\vee(1-q_{k})>0.$$ (A.12)

The last term in equation (A.3.1) is greater than

$$p_{k}\min_{m}n_{k,m}\big{(}E[A^{2}_{k,m}(1-A_{k,m})^{2}]-(E[A_{k,m}(1-A_{k,m})])^{2}\big{)}\sum_{m=1}^{m_{k}}\frac{n_{k,m}}{n_{k}}(\tau_{k,m}-\tau_{k})^{2},$$

which converges to infinity because $p_{k}\min_{m}n_{k,m}\rightarrow\infty$. That is, the last term dominates the variance in large samples provided that (A.11) and (A.12) hold.

We will now derive the large sample distribution of $\widehat{\tau}_{k}^{\rm{\,fixed}}$. To show that Lyapunov’s condition holds for $F_{k}$, notice that

$$|( \psi_{k,m}-\hbox{\vbox{\hrule height=0.5pt\kern 0.77498pt\hbox{\kern-0.18004pt$\psi$\kern-0.18004pt}}}_{k,m})+(\varphi_{k,m}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\varphi$\kern-0.18004pt}}}_{k,m})|^{3}$$

$$=\frac{1}{(n_{k}p_{k}q_{k})^{3/2}}\Bigg{|}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}((e_{k,i}(1)+\tau_{k,m}-\tau_{k})W_{k,i}+e_{k,i}(0)(1-W_{k,i}))(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})$$

$$\qquad\qquad\qquad\qquad\qquad-(\tau_{k,m}-\tau)q_{k}E[A_{k,m}(1-A_{k,m})](1-(1-p_{k})^{n_{k,m}})\Bigg{|}^{3},$$

where the last term inside the absolute value comes from the bias correction. Notice that,

$$\Bigg{|}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\} R_{k,i}(e_{k,i}(1)+\tau_{k,m}-\tau_{k})W_{k,i}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})\Bigg{|}^{3}$$

$$=\Bigg{|}(1-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}(e_{k,i}(1)+\tau_{k,m}-\tau_{k})W_{k,i}\Bigg{|}^{3}$$

$$\leq c\,\Bigg{|}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}\Bigg{|}^{3}$$

$$\leq c\,\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$N$\kern-0.18004pt}}}_{k,m}^{3}.$$

From the formula of the third moment of a binomial random variable, we obtain

$$E[\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$N$\kern-0.18004pt}}}_{k,m}^{3}] =q_{k}E[\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$N$\kern-0.18004pt}}}_{k,m}^{3}|Q_{k,m}=1]$$

$$=n_{k,m}^{3}p_{k}^{3}q_{k}+\scaleto{\mathcal{O}}{5pt}(n_{k,m}^{3}p_{k}^{3}q_{k}),$$

as $p_{k}n_{k,m}\rightarrow\infty$. Now,

$$\frac{1}{f_{k}^{3/2}}\sum_{k=1}^{m_{k}}E\Bigg{[}\Bigg{|} \frac{1}{\sqrt{n_{k}p_{k}q_{k}}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}(e_{k,i}(1)+\tau_{k,m}-\tau_{k})W_{k,i}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})\Bigg{|}^{3}\Bigg{]}$$

$$\leq c\,\frac{n_{k}\max_{n}n_{k,m}^{2}p^{3}_{k}q_{k}}{(n_{k}p_{k}q_{k})^{3/2}(p_{k}\min_{m}n_{k,m})^{3/2}}=c\left(\frac{\max_{m}n_{k,m}}{\min_{m}n_{k,m}}\right)^{2}\frac{1}{(m_{k}q_{k})^{1/2}}\longrightarrow 0.$$

Similar calculations deliver the analogous result for the term involving $e_{k,i}(0)$, and proving the result for the bias term is straightforward. Therefore, we obtain

$$\frac{1}{f_{k}^{3/2}}\sum_{m=1}^{m_{k}}|(\psi_{k,m}-\hbox{\vbox{\hrule height=0.5pt\kern 0.77498pt\hbox{\kern-0.18004pt$\psi$\kern-0.18004pt}}}_{k,m})+(\varphi_{k,m}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\varphi$\kern-0.18004pt}}}_{k,m})|^{3}\longrightarrow 0.$$

By the Central Limit Theorem for arrays, this implies

$$\sqrt{n_{k}p_{k}q_{k}}F_{k}/f_{k}^{1/2}\stackrel{{\scriptstyle d}}{{\longrightarrow}}N(0,1).$$

Let $\tilde{v}_{k}=f_{k}/(\mu_{k}(1-\mu_{k})-\sigma^{2}_{k})^{2}$. Then,

$$\sqrt{n_{k}p_{k}q_{k}}(\widehat{\tau}_{k}^{\rm{\,fixed}}-\tau_{k})/\tilde{v}_{k}^{1/2}\stackrel{{\scriptstyle d}}{{\longrightarrow}}N(0,1).$$

As a result,

$$\sqrt{N_{k}}(\widehat{\tau}_{k}^{\rm{\,fixed}}-\tau_{k})/\tilde{v}_{k}^{1/2}\stackrel{{\scriptstyle d}}{{\longrightarrow}}N(0,1).$$

A.3.2 Estimation of the variance

Let

$$N_{k,m,0} =\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}(1-W_{k,i})$$

and

$$N_{k,m,1} =\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}.$$

Let

$$\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$Y$\kern-0.18004pt}}}_{k,m}=\frac{1}{\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$N$\kern-0.18004pt}}}_{k,m}\vee 1}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}Y_{k,i}.$$

Then,

$$\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$Y$\kern-0.18004pt}}}_{k,m}=\widehat{\alpha}_{k,m}+\widehat{\tau}_{k,m}\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m},$$

where

$$\widehat{\alpha}_{k,m}=\frac{1}{N_{k,m,0}\vee 1}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}(1-W_{k,i})Y_{k,i},$$

$$\widehat{\tau}_{k,m}=\frac{1}{N_{k,m,1}\vee 1}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}Y_{k,i}-\frac{1}{N_{k,m,0}\vee 1}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}(1-W_{k,i})Y_{k,i},$$

and, as before,

$$\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}=\frac{1}{\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$N$\kern-0.18004pt}}}_{k,m}\vee 1}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}.$$

Let $\widetilde{U}_{k,i}=\widetilde{Y}_{k,i}-\widehat{\tau}_{k}^{\,{\rm fixed}}\widetilde{W}_{k,i}$, where $\widetilde{Y}_{k,i}=Y_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$Y$\kern-0.18004pt}}}_{k,m_{k,i}}$, $\widetilde{W}_{k,i}=(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m_{k,i}})$, and $\widehat{\tau}_{k}^{\,\rm fixed}$ is the within estimator of $\tau_{k}$. Let $\widetilde{\Sigma}_{k}=\sum_{m=1}^{m_{k}}\widetilde{\Sigma}_{k,m}$, where

$$\widetilde{\Sigma}_{k,m} =\left(\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}\widetilde{W}_{k,i}\widetilde{U}_{k,i}\right)^{2}.$$

Also, let

$$\widetilde{Q}_{k}=\sum_{i=1}^{n_{k}}R_{k,i}\widetilde{W}_{k,i}^{2}.$$

Then, the cluster estimator of the variance of $\sqrt{N_{k}}(\widehat{\tau}^{\,{\rm fixed}}_{k}-\tau_{k})$ is

$$\widetilde{V}_{k}^{\rm{cluster}}=N_{k}\widetilde{Q}_{k}^{-1}\widetilde{\Sigma}_{k}\widetilde{Q}_{k}^{-1}.$$

We know already that

$$\frac{1}{n_{k}p_{k}q_{k}}\widetilde{Q}_{k}-(\mu_{k}(1-\mu_{k})-\sigma_{k}^{2})\stackrel{{\scriptstyle p}}{{\longrightarrow}}0,$$

with $\mu_{k}(1-\mu_{k})-\sigma_{k}^{2}$ bounded away from zero. To establish convergence of $\widetilde{\Sigma}_{k}/(n_{k}p_{k}q_{k}f_{k})$, first notice that, for $m_{k,i}=m$, we have

$$\widetilde{U}_{k,i} =Y_{k,i}-(\widehat{\alpha}_{k,m}+\widehat{\tau}_{k,m}\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})-\widehat{\tau}_{k}^{\,{\rm fixed}}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})$$

$$=y_{k,i}(1)W_{k,i}+y_{k,i}(0)(1-W_{k,i})-(\alpha_{k,m}+\tau_{k,m}\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})-\widehat{\tau}_{k}^{\,{\rm fixed}}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})$$

$$-(\widehat{\alpha}_{k,m}-\alpha_{k,m})-(\widehat{\tau}_{k,m}-\tau_{k,m})\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}$$

$$=e_{k,i}(1)W_{k,i}+e_{k,i}(0)(1-W_{k,i})+(\tau_{k,m}-\widehat{\tau}_{k}^{\,{\rm fixed}})(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})$$

$$-(\widehat{\alpha}_{k,m}-\alpha_{k,m})-(\widehat{\tau}_{k,m}-\tau_{k,m})\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}$$

$$=e_{k,i}(1)W_{k,i}+e_{k,i}(0)(1-W_{k,i})+(\tau_{k,m}-\tau_{k})(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})$$

$$-(\widehat{\tau}_{k}^{\,{\rm fixed}}-\tau_{k})(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})-(\widehat{\alpha}_{k,m}-\alpha_{k,m})-(\widehat{\tau}_{k,m}-\tau_{k,m})\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}.$$

For $m_{k,i}=m$ and $N_{k,m,0},N_{k,m,1}\geq 1$, let

$$\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$U$\kern-0.18004pt}}}_{k,i}=e_{k,i}(1)W_{k,i}+e_{k,i}(0)(1-W_{k,i})+(\tau_{k,m}-\tau_{k})(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}),$$

and let $\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$U$\kern-0.18004pt}}}_{k,i}=0$ for $m_{k,i}=m$ and $N_{k,m,0}N_{k,m,1}=0$. Then, for $m_{k,i}=m$ and $N_{k,m,0}N_{k,m,1}\geq 1$, we have

$$\widetilde{U}_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$U$\kern-0.18004pt}}}_{k,i}=-(\widehat{\tau}_{k}^{\,{\rm fixed}}-\tau_{k})(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})-(\widehat{\alpha}_{k,m}-\alpha_{k,m})-(\widehat{\tau}_{k,m}-\tau_{k,m})\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}.$$

Then,

$$\Bigg{(} \sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}\widetilde{W}_{k,i}\widetilde{U}_{k,i}\Bigg{)}^{2}$$

$$=\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}\widetilde{W}_{k,i}\Big{(}\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$U$\kern-0.18004pt}}}_{k,i}+\big{(}\widetilde{U}_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$U$\kern-0.18004pt}}}_{k,i}\big{)}\Big{)}\Bigg{)}^{2}$$

$$=\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}\widetilde{W}_{k,i}\Big{(}\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$U$\kern-0.18004pt}}}_{k,i}-(\widehat{\tau}_{k}^{\,{\rm fixed}}-\tau_{k})(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})\Big{)}\Bigg{)}^{2}$$

$$=\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}\widetilde{W}_{k,i}\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$U$\kern-0.18004pt}}}_{k,i}-(\widehat{\tau}_{k}^{\,{\rm fixed}}-\tau_{k})\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})\Bigg{)}^{2}.$$

Using the formula for the second moment of a binomial distribution and $n\geq 1$, we obtain,

$$E\Bigg{[}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i} W_{k,i}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})\Bigg{)}^{2}\Big{|}\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$N$\kern-0.18004pt}}}_{k,m}=n\Bigg{]}$$

$$=E\Bigg{[}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}(1-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})R_{k,i}W_{k,i}\Bigg{)}^{2}\Big{|}\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$N$\kern-0.18004pt}}}_{k,m}=n\Bigg{]}$$

$$\leq E\Bigg{[}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}\Bigg{)}^{2}\Big{|}\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$N$\kern-0.18004pt}}}_{k,m}=n\Bigg{]}$$

$$\leq n^{2}+n.$$

From the formula of the sum of the first two moments of a binomial distribution, we obtain

$$\sum_{m=1}^{m_{k}}E\Bigg{[}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})\Bigg{)}^{2}\Bigg{]}\leq\sum_{m=1}^{m_{k}}(n_{k,m}^{2}p_{k}^{2}q_{k}+n_{k,m}p_{k}q_{k}).$$

Therefore,

$$\frac{1}{n_{k}p_{k}q_{k}f_{k}}(\widehat{\tau}_{k}^{\,{\rm fixed}}-\tau_{k})^{2} \sum_{m=1}^{m_{k}}E\Bigg{[}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})\Bigg{)}^{2}\Bigg{]}$$

$$\leq\frac{n_{k}p_{k}q_{k}}{f_{k}}(\widehat{\tau}_{k}^{\,{\rm fixed}}-\tau_{k})^{2}\frac{1}{(n_{k}p_{k}q_{k})^{2}}\sum_{m=1}^{m_{k}}(n_{k,m}^{2}p_{k}^{2}q_{k}+n_{k,m}p_{k}q_{k})$$

$$=\mathcal{O}_{p}(1)\left(\frac{\max_{m}n_{k,m}}{\min_{m}n_{k,m}}\frac{1}{m_{k}q_{k}}+\frac{1}{n_{k}p_{k}q_{k}}\right)\stackrel{{\scriptstyle p}}{{\longrightarrow}}0.$$

Now, notice that

$$\frac{1}{n_{k}p_{k}q_{k}f_{k}} \sum_{m=1}^{m_{k}}\left(\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}\widetilde{W}_{k,i}\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$U$\kern-0.18004pt}}}_{k,i}\right)^{2}$$

$$=\frac{1}{n_{k}p_{k}q_{k}f_{k}}\sum_{m=1}^{m_{k}}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}\big{(}e_{k,i}(1)W_{k,i}+e_{k,i}(0)(1-W_{k,i})\big{)}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})$$

$$\hskip 170.71652pt+(\tau_{k,m}-\tau_{k})\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})^{2}\Bigg{)}^{2}.$$

Equation (A.9) (and the analogous result for the sum involving terms with $e_{k,i}(0)$), implies

$$\frac{1}{n_{k}p_{k}q_{k}f_{k}}\sum_{m=1}^{m_{k}}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}\big{(}e_{k,i}(1)W_{k,i}+e_{k,i}(0)(1-W_{k,i})\big{)}(\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}-A_{k,m})\Bigg{)}^{2}\stackrel{{\scriptstyle p}}{{\longrightarrow}}0.$$

As a result, it is enough to establish convergence of $\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\Sigma$\kern-0.18004pt}}}_{k}/(n_{k}p_{k}q_{k}f_{k})$, where

$$\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\Sigma$\kern-0.18004pt}}}_{k} =\sum_{m=1}^{m_{k}}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}\big{(}e_{k,i}(1)W_{k,i}+e_{k,i}(0)(1-W_{k,i})\big{)}(W_{k,i}-A_{k,m})$$

$$+(\tau_{k,m}-\tau_{k})\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})^{2}\Bigg{)}^{2}$$

$$=\sum_{m=1}^{m_{k}}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}\Big{(}R_{k,i}W_{k,i}(W_{k,i}-A_{k,m})-p_{k}q_{k}A_{k,m}(1-A_{k,m})\Big{)}e_{k,i}(1)$$

$$+\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}\Big{(}R_{k,i}(1-W_{k,i})(W_{k,i}-A_{k,m})+p_{k}q_{k}A_{k,m}(1-A_{k,m})\Big{)}e_{k,i}(0)$$

$$+(\tau_{k,m}-\tau_{k})\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})^{2}\Bigg{)}^{2}.$$

We will next show that

$$\frac{1}{n_{k}p_{k}q_{k}f_{k}}\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\Sigma$\kern-0.18004pt}}}_{k}-\frac{f_{k}^{\rm cluster}}{f_{k}}\stackrel{{\scriptstyle p}}{{\longrightarrow}}0,$$ (A.13)

where

$$f_{k}^{\rm cluster} =\frac{1}{n_{k}}E[A_{k,m}(1-A_{k,m})^{2}]\sum_{i=1}^{n_{k}}e_{k,i}^{2}(1)$$

$$+\frac{1}{n_{k}}E[A^{2}_{k,m}(1-A_{k,m})]\sum_{i=1}^{n_{k}}e_{k,i}^{2}(0)$$

$$-\frac{1}{n_{k}}p_{k}E[A^{2}_{k,m}(1-A_{k,m})^{2}]\sum_{i=1}^{n_{k}}(e_{k,i}(1)-e_{k,i}(0))^{2}$$

$$+(E[A_{k,m}(1-A_{k,m})]-(5+p_{k})E[A^{2}_{k,m}(1-A_{k,m})^{2}])\sum_{m=1}^{m_{k}}\frac{n_{k,m}}{n_{k}}(\tau_{k,m}-\tau_{k})^{2}$$

$$+p_{k}E[A^{2}_{k,m}(1-A_{k,m})^{2}]\sum_{m=1}^{m_{k}}\frac{n^{2}_{k,m}}{n_{k}}(\tau_{k,m}-\tau_{k})^{2}.$$

Let

$$X_{k,m} =\frac{1}{n_{k}p_{k}q_{k}}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}\big{(}e_{k,i}(0)(1-W_{k,i})$$

$$+e_{k,i}(1)W_{k,i}\big{)}(W_{k,i}-A_{k,m})+(\tau_{k,m}-\tau_{k})\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})^{2}\Bigg{)}^{2}$$

Using the result in equation (A.3.1) and results on the moments of the binomial distribution (see intermediate calculations in section A.7), we obtain

$$\frac{1}{n_{k}p_{k}q_{k}}E[\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$\Sigma$\kern-0.18004pt}}}_{k}] =\sum_{m=1}^{m_{k}}E[X_{k,m}]$$

$$=f_{k}^{\rm cluster}+\scaleto{\mathcal{O}}{5pt}(1).$$

Therefore, to show that equation (A.13) holds, we will show

$$\frac{1}{f_{k}^{2}}\sum_{m=1}^{m_{k}}E[X_{k,m}^{2}]\longrightarrow 0.$$ (A.14)

Let

$$\theta_{k} =E[(R_{k,i}W_{k,i}(W_{k,i}-A_{k,m})-p_{k}A_{k,m}(1-A_{k,m}))^{2}|m_{k,i}=m,Q_{k,m}=1]$$

$$=p_{k}\left(E[A_{k,m}(1-A_{k,m})^{2}]-p_{k}E[A_{k,m}^{2}(1-A_{k,m})^{2}]\right),$$

and

$$\pi_{k} =E[(R_{k,i}W_{k,i}(W_{k,i}-A_{k,m})-p_{k}A_{k,m}(1-A_{k,m}))^{4}|m_{k,i}=m,Q_{k,m}=1]$$

$$=p_{k}E[(W_{k,i}(W_{k,i}-A_{k,m})-p_{k}A_{k,m}(1-A_{k,m}))^{4}|m_{k,i}=m]+p_{k}^{4}(1-p_{k})E[A^{4}_{k,m}(1-A_{k,m})^{4}].$$

Let

$$X_{k,m,1} =\frac{1}{n_{k}p_{k}q_{k}}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}(W_{k,i}-A_{k,m})e_{k,i}(1)\Bigg{)}^{2}$$

$$=\frac{Q_{k,m}}{n_{k}p_{k}q_{k}}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}(R_{k,i}W_{k,i}(W_{k,i}-A_{k,m})-p_{k}A_{k,m}(1-A_{k,m}))e_{k,i}(1)\Bigg{)}^{2}.$$

Then,

$$E[X_{k,m,1}^{2}] =q_{k}E[X_{k,m,1}^{2}|Q_{k,m}=1]$$

$$=\frac{\pi_{k}}{n^{2}_{k}p^{2}_{k}q_{k}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}e_{k,i}^{4}(1)$$

$$+\frac{6\theta^{2}_{k}}{n^{2}_{k}p^{2}_{k}q_{k}}\sum_{i=1}^{n_{k}-1}\sum_{j=i+1}^{n_{k}}1\{m_{k,i}=m_{k,j}=m\}e_{k,i}^{2}(1)e_{k,j}^{2}(1).$$

Therefore, because $n_{k}p_{k}q_{k}\rightarrow\infty$ and $m_{k}q_{k}\rightarrow\infty$, we obtain

$$\sum_{m=1}^{m_{k}}E[X_{k,m,1}^{2}] \leq\frac{c}{n_{k}p_{k}q_{k}}\left(\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}e_{k,i}^{4}(1)\right)+\frac{c}{m_{k}q_{k}}\frac{\max_{m}n^{2}_{k,m}}{\min_{m}n_{k,m}^{2}}$$

$$\hskip 28.45274pt\times\left(\frac{1}{m_{k}}\sum_{m=1}^{m_{k}}\frac{1}{\max_{m}n^{2}_{k,m}}\sum_{i=1}^{n_{k}-1}\sum_{j=i+1}^{n_{k}}1\{m_{k,i}=m_{k,j}=m\}e_{k,i}^{2}(1)e_{k,j}^{2}(1)\right)$$

$$\longrightarrow 0.$$ (A.15)

Using the same argument, we obtain

$$\sum_{m=1}^{m_{k}}E[X_{k,m,2}^{2}]\longrightarrow 0,$$ (A.16)

where

$$X_{k,m,2}=\frac{1}{n_{k}p_{k}q_{k}}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}(1-W_{k,i})(W_{k,i}-A_{k,m})e_{k,i}(0)\Bigg{)}^{2}.$$

Notice that equations (A.3.2) and (A.16) imply

$$\frac{1}{f_{k}^{2}}\sum_{m=1}^{m_{k}}E[X_{k,m,1}^{2}]\longrightarrow 0$$

and

$$\frac{1}{f_{k}^{2}}\sum_{m=1}^{m_{k}}E[X_{k,m,2}^{2}]\longrightarrow 0.$$

Notice that the last two equations hold even if $f_{k}$ is bounded (e.g., when $\tau_{k,m}-\tau_{k}=0$ for all $k$ and $m$), as long as $f_{k}$ is bounded away from zero in large samples. In section A.3.3 we derive conditions so that $f_{k}$ is bounded away from zero in large samples even if $\tau_{k,m}-\tau_{k}=0$ for all $k$ and $m$. Now, let

$$X_{k,m,3}=\frac{1}{n_{k}p_{k}q_{k}}\Bigg{(}(\tau_{k,m}-\tau_{k})\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})\Bigg{)}^{2}.$$

Recall that, under the conditions in (A.11) and (A.12), $f_{k}\rightarrow\infty$ and $p_{k}\min n_{k,m}/f_{k}$ is bounded for large $k$ and, therefore, $p_{k}\max n_{k,m}/f_{k}$ is bounded for large $k$. Then (see intermediate calculations at the end of this document), for large $k$,

$$\frac{1}{f_{k}^{2}}\sum_{m=1}^{m_{k}}E[X_{k,m,3}^{2}] =\frac{1}{(n_{k}p_{k}q_{k}f_{k})^{2}}\sum_{m=1}^{m_{k}}n_{k,m}^{4}p_{k}^{4}q_{k}(\tau_{k,m}-\tau_{k})^{4}\left(1+\mathcal{O}\left(\frac{1}{p_{k}\min_{m}n_{k,m}}\right)\right)$$

$$=\frac{p_{k}\max_{m}n^{2}_{k,m}}{m_{k}q_{k}f_{k}\min_{m}n_{k,m}}\left(\frac{p_{k}}{f_{k}}\sum_{m=1}^{m_{k}}\frac{n_{k,m}^{2}}{n_{k}}(\tau_{k,m}-\tau_{k})^{4}\right)\left(1+\mathcal{O}\left(\frac{1}{p_{k}\min_{m}n_{k,m}}\right)\right)$$

$$=\mathcal{O}\left(\frac{1}{m_{k}q_{k}}\right)\left(1+\mathcal{O}\left(\frac{1}{p_{k}\min_{m}n_{k,m}}\right)\right)\rightarrow 0.$$

Now, Hölder’s inequality implies that equation (A.14) holds (see intermediate calculations).

Now let,

$$\tilde{v}_{k}^{\rm cluster}=f_{k}^{\rm cluster}/(\mu_{k}(1-\mu_{k})-\sigma^{2}_{k})^{2}.$$

We obtain,

$$\frac{\widetilde{V}_{k}^{\rm{cluster}}}{\tilde{v}_{k}}=\frac{\tilde{v}_{k}^{\rm cluster}}{\tilde{v}_{k}}+\scaleto{\mathcal{O}}{5pt}_{p}(1).$$

We will next establish the analogous result for the heteroskedaticity-robust variance estimator. Let

$$\widetilde{\Sigma}_{k}^{\rm robust} =\sum_{i=1}^{n_{k}}R_{k,i}\widetilde{W}_{k,i}^{2}\widetilde{U}_{k,i}^{2}.$$

Then, the heteroskedasticity-robust estimator of the variance of $\sqrt{N_{k}}(\widehat{\tau}^{\,{\rm fixed}}_{k}-\tau_{k})$ is

$$\widetilde{V}_{k}^{\rm{robust}}=N_{k}\widetilde{Q}_{k}^{-1}\widetilde{\Sigma}_{k}^{\rm robust}\widetilde{Q}_{k}^{-1}.$$

As we have established before,

$$\widetilde{U}_{k,i} =e_{k,i}(1)W_{k,i}+e_{k,i}(0)(1-W_{k,i})+(\tau_{k,m}-\tau_{k})(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})$$

$$-(\widehat{\tau}_{k}^{\,{\rm fixed}}-\tau_{k})(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})-(\widehat{\alpha}_{k,m}-\alpha_{k,m})-(\widehat{\tau}_{k,m}-\tau_{k,m})\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}.$$

For $m_{k,i}=m$ and $N_{k,m,0}N_{k,m,1}\geq 1$, let

$$\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$U$\kern-0.18004pt}}}_{k,i}=e_{k,i}(1)W_{k,i}+e_{k,i}(0)(1-W_{k,i})+(\tau_{k,m}-\tau_{k})(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}),$$

and let $\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$U$\kern-0.18004pt}}}_{k,i}=0$ for $m_{k,i}=m$ and $N_{k,m,0}N_{k,m,1}=0$. Then, for $m_{k,i}=m$ and $N_{k,m,0}N_{k,m,1}\geq 1$, we have

$$\widetilde{U}_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$U$\kern-0.18004pt}}}_{k,i}=-(\widehat{\tau}_{k}^{\,{\rm fixed}}-\tau_{k})(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})-(\widehat{\alpha}_{k,m}-\alpha_{k,m})-(\widehat{\tau}_{k,m}-\tau_{k,m})\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m},$$

and

$$\frac{1}{n_{k}p_{k}q_{k}}\sum_{i=1}^{n_{k}}R_{k,i}\widetilde{W}_{k,i}^{2}\widetilde{U}_{k,i}^{2}=\frac{1}{n_{k}p_{k}q_{k}}\sum_{i=1}^{n_{k}}R_{k,i}\widetilde{W}_{k,i}^{2}\Big{(}\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$U$\kern-0.18004pt}}}_{k,i}+\big{(}\widetilde{U}_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$U$\kern-0.18004pt}}}_{k,i}\big{)}\Big{)}^{2}.$$ (A.17)

Focusing on the part of the right hand side of last equation that depends on the first term of $\widetilde{U}_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$U$\kern-0.18004pt}}}_{k,i}$, we obtain

$$\frac{1}{n_{k}p_{k}q_{k}}\sum_{i=1}^{n_{k}}R_{k,i}\widetilde{W}_{k,i}^{4}(\widehat{\tau}_{k}^{\,{\rm fixed}}-\tau_{k})^{2}\leq(\widehat{\tau}_{k}^{\,{\rm fixed}}-\tau_{k})^{2}\frac{1}{n_{k}p_{k}q_{k}}\sum_{i=1}^{n_{k}}R_{k,i}\widetilde{W}_{k,i}^{2}\stackrel{{\scriptstyle p}}{{\longrightarrow}}0.$$

We will focus now on the part of the right-hand side of equation (A.17) that that depends on the second term of $\widetilde{U}_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$U$\kern-0.18004pt}}}_{k,i}$,

$$\frac{1}{n_{k}p_{k}q_{k}} \sum_{m=1}^{m_{k}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}\widetilde{W}_{k,i}^{2}(\widehat{\alpha}_{k,m}-\alpha_{k,m})^{2}.$$

Using the formula for the variance of a sample mean under sampling without replacement [e.g., in the supplement of abadie2020sampling], we obtain for $1\leq n\leq n_{k,m}-1$,

$$E\Bigg{[}(\widehat{\alpha}_{k,m}-\alpha_{k,m})^{2}\sum_{i=1}^{n_{k}} 1\{m_{k,i}=m\}R_{k,i}\widetilde{W}_{k,i}^{2}\Big{|}N_{k,m,0}=n\Bigg{]}$$

$$=E\Bigg{[}(\widehat{\alpha}_{k,m}-\alpha_{k,m})^{2}\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$N$\kern-0.18004pt}}}_{k,m}\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}(1-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})\Big{|}N_{k,m,0}=n\Bigg{]}$$

$$\leq E\Big{[}n(\widehat{\alpha}_{k,m}-\alpha_{k,m})^{2}\big{|}N_{k,m,0}=n\Big{]}$$

$$=n\,\mbox{var}(\widehat{\alpha}_{k,m}|N_{k,m,0}=n)$$

$$=s^{2}_{k,m,0}\Big{(}1-\frac{n}{n_{k,m}}\Big{)},$$ (A.18)

where

$$s^{2}_{k,m,0}=\frac{1}{n_{k,m}-1}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}(y_{k,i}(0)-\alpha_{k,m})^{2}.$$

Because $s^{2}_{k,m,0}$ is bounded, so is the right-hand side of equation (A.3.2). As a result

$$E\Bigg{[}\frac{1}{n_{k}p_{k}q_{k}}\sum_{m=1}^{m_{k}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}\widetilde{W}_{k,i}^{2}(\widehat{\alpha}_{k,m}-\alpha_{k,m})^{2}\Bigg{]}\leq c\,\frac{m_{k}}{n_{k}p_{k}}\longrightarrow 0.$$

An analogous derivation applies to the part of the right-hand side of equation (A.17) that depends on the third term of $\widetilde{U}_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$U$\kern-0.18004pt}}}_{k,i}$. (Notice that $\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}\leq 1$ and that $\widehat{\tau}_{k,m}-\tau_{k,m}$ is equal to minus the difference between $\widehat{\alpha}_{k,m}-\alpha_{k,m}$ and the analogous difference for the treated.

Therefore, we will study the behavior of

$$\frac{1}{n_{k}p_{k}q_{k}}\sum_{i=1}^{n_{k}}R_{k,i}\widetilde{W}_{k,i}^{2}\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$U$\kern-0.18004pt}}}_{k,i}^{2}.$$ (A.19)

First, notice that

$$\frac{1}{n_{k}p_{k}q_{k}} \sum_{m=1}^{m_{k}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}\Big{|}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})^{2}-(W_{k,i}-A_{k,m})^{2}\Big{|}W_{k,i}e^{2}_{k,i}(1)$$

$$\leq c\frac{1}{n_{k}p_{k}q_{k}}\sum_{m=1}^{m_{k}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}\Big{|}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})+(W_{k,i}-A_{k,m})\Big{|}\big{|}\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}-A_{k,m}\big{|}W_{k,i}$$

$$\leq c\Bigg{(}\frac{1}{n_{k}p_{k}q_{k}}\sum_{m=1}^{m_{k}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}\Big{(}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})+(W_{k,i}-A_{k,m})\Big{)}^{2}W_{k,i}\Bigg{)}^{1/2}$$

$$\quad\times\Bigg{(}\frac{1}{n_{k}p_{k}q_{k}}\sum_{m=1}^{m_{k}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}\big{(}\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}-A_{k,m}\big{)}^{2}\Bigg{)}^{1/2}.$$ (A.20)

The inside of the first square root in equation (A.3.2) is bounded by a constant times

$$\frac{1}{n_{k}p_{k}q_{k}}\sum_{m=1}^{m_{k}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i},$$

which converges in probability to one. The expectation of the inside of the second square root in equation (A.3.2) is

$$\frac{1}{n_{k}p_{k}q_{k}}\sum_{m=1}^{m_{k}}E\big{[}\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$N$\kern-0.18004pt}}}_{k,m}\big{(}\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}-A_{k,m}\big{)}^{2}\big{]}\leq c\frac{m_{k}}{n_{k}p_{k}}\longrightarrow 0.$$

As a result, the right-hand side of equation (A.3.2) converges to zero in probability. The derivation with $(1-W_{k,i})e^{2}_{k,i}(0)$ replacing $W_{k,i}e^{2}_{k,i}(1)$ in equation (A.3.2) is analogous. Now, notice that

$$(W_{k,i} -\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})^{4}-(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$A$\kern-0.18004pt}}}_{k,m})^{4}$$

$$=-\big{(}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})^{2}+(W_{k,i}-A_{k,m})^{2}\big{)}\big{(}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})+(W_{k,i}-A_{k,m})\big{)}(\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}-A_{k,m}).$$

Because the first factor of the expression above is bounded, we obtain

$$\frac{1}{n_{k}p_{k}q_{k}} \sum_{m=1}^{m_{k}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}\Big{|}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})^{4}-(W_{k,i}-A_{k,m})^{4}\Big{|}(\tau_{k,m}-\tau_{k})^{2}$$

$$\leq c\Bigg{(}\frac{1}{n_{k}p_{k}q_{k}}\sum_{m=1}^{m_{k}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}\Bigg{)}^{1/2}$$

$$\quad\times\Bigg{(}\frac{1}{n_{k}p_{k}q_{k}}\sum_{m=1}^{m_{k}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}\big{(}\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}-A_{k,m}\big{)}^{2}\Bigg{)}^{1/2}.$$ (A.21)

Now, the right-hand side of equation (A.3.2) converges to zero in probability by the same argument as for equation (A.3.2). Cauchy-Schwarz inequality implies,

$$\frac{1}{n_{k}p_{k}q_{k}}\sum_{i=1}^{n_{k}}R_{k,i}\widetilde{W}_{k,i}^{2}\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$U$\kern-0.18004pt}}}_{k,i}^{2}=\frac{1}{n_{k}p_{k}q_{k}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}(W_{k,i}-A_{k,m})^{2}\breve{U}_{k,i}^{2}+\scaleto{\mathcal{O}}{5pt}_{p}(1),$$

where

$$\breve{U}_{k,i}=e_{k,i}(1)W_{k,i}+e_{k,i}(0)(1-W_{k,i})+(\tau_{k,m}-\tau_{k})(W_{k,i}-A_{k,m}),$$ (A.22)

for $m_{k,i}=m$ and $N_{k,m,0}N_{k,m,1}\geq 1$, and $\breve{U}_{k,i}=0$ for $N_{k,m,0}N_{k,m,1}\geq 0$. Therefore, we will study the behavior of

$$\frac{1}{n_{k}p_{k}q_{k}}\sum_{m=1}^{m_{k}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}(W_{k,i}-A_{k,m})^{2}\breve{U}_{k,i}^{2}.$$

We know,

$$\frac{1}{n_{k}p_{k}q_{k}}\sum_{m=1}^{m_{k}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i} (W_{k,i}-A_{k,m})^{2}W_{k,i}e^{2}_{k,i}(1)$$

$$-E[A_{k,m}(1-A_{k,m})^{2}]\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}e^{2}_{k,i}(1)\stackrel{{\scriptstyle p}}{{\longrightarrow}}0,$$

and

$$\frac{1}{n_{k}p_{k}q_{k}}\sum_{m=1}^{m_{k}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i} (W_{k,i}-A_{k,m})^{2}(1-W_{k,i})e^{2}_{k,i}(0)$$

$$-E[A^{2}_{k,m}(1-A_{k,m})]\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}e^{2}_{k,i}(0)\stackrel{{\scriptstyle p}}{{\longrightarrow}}0.$$

Now, notice that

$$E[(W_{k,i}-A_{k,m})^{4}|m_{k,i}=m,R_{k,i}=1,A_{k,m}=a] =(1-a)^{4}a+a^{4}(1-a)$$

$$=a(1-a)[(1-a)^{3}+a^{3}]$$

$$=a(1-a)[1-3a(1-a)],$$

which implies

$$E\Bigg{[}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i} (W_{k,i}-A_{k,m})^{4}(\tau_{k,m}-\tau_{k})^{2}\Bigg{]}$$

$$=n_{k,m}p_{k}q_{k}E[A_{k,m}(1-A_{k,m})(1-3A_{k,m}(1-A_{k,m}))](\tau_{k,m}-\tau_{k})^{2},$$

and

$$E\Bigg{[}\frac{1}{n_{k}p_{k}q_{k}}\sum_{m=1}^{m_{k}} \sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}(W_{k,i}-A_{k,m})^{4}(\tau_{k,m}-\tau_{k})^{2}\Bigg{]}$$

$$=E[A_{k,m}(1-A_{k,m})(1-3A_{k,m}(1-A_{k,m}))]\sum_{m=1}^{m_{k}}\frac{n_{k,m}}{n_{k}}(\tau_{k,m}-\tau_{k})^{2}.$$

Notice now that

$$\frac{1}{(n_{k}p_{k}q_{k})^{2}}\sum_{m=1}^{m_{k}} E\Bigg{[}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}(W_{k,i}-A_{k,m})^{4}(\tau_{k,m}-\tau_{k})^{2}\Bigg{)}^{2}\Bigg{]}$$

$$\leq c\,\frac{1}{(n_{k}p_{k}q_{k})^{2}}\sum_{m=1}^{m_{k}}E\Bigg{[}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}\Bigg{)}^{2}\Bigg{]}$$

$$\leq c\,\frac{q_{k}}{(n_{k}p_{k}q_{k})^{2}}\sum_{m=1}^{m_{k}}(n_{k,m}p_{k}+n_{k,m}^{2}p_{k}^{2})$$

$$=c\Bigg{(}\frac{1}{n_{k}p_{k}q_{k}}+\frac{\max_{m}n_{k,m}}{\min_{m}n_{k,m}}\frac{1}{m_{k}q_{k}}\Bigg{)}\stackrel{{\scriptstyle p}}{{\longrightarrow}}0.$$

Notice also that expectations of the sums of products of the terms on the right-hand side of equation (A.22) are equal to zero. Then,

$$\frac{1}{n_{k}p_{k}q_{k}}\widetilde{\Sigma}_{k}^{\rm robust}-f_{k}^{\rm robust}\stackrel{{\scriptstyle p}}{{\longrightarrow}}0,$$

where

$$f_{k}^{\rm robust} =E[A_{k,m}(1-A_{k,m})^{2}]\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}e^{2}_{k,i}(1)+E[A^{2}_{k,m}(1-A_{k,m})]\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}e^{2}_{k,i}(0)$$

$$+E[A_{k,m}(1-A_{k,m})(1-3A_{k,m}(1-A_{k,m}))]\sum_{m=1}^{m_{k}}\frac{n_{k,m}}{n_{k}}(\tau_{k,m}-\tau_{k})^{2}.$$

Now let,

$$\tilde{v}_{k}^{\rm robust}=f_{k}^{\rm robust}/(\mu_{k}(1-\mu_{k})-\sigma^{2}_{k})^{2}.$$

We obtain,

$$\widetilde{V}_{k}^{\rm{robust}}=\tilde{v}_{k}^{\rm robust}+\scaleto{\mathcal{O}}{5pt}_{p}(1).$$

A.3.3 Large $k$ results the fixed effects case under homogeneous average treatment effects across clusters

We will now study the Lyapounov’s condition for the case $\tau_{k,m}=\tau_{k}$ for all $k$ and $m=1,\ldots,m_{k}$, so

$$f_{k}=\sum_{m=1}^{m_{k}}E[\psi_{k,m}^{2}].$$

Notice that

$$\sum_{m=1}^{m_{k}}E[\psi_{k,m}^{2}] \geq\frac{1}{n_{k}}E[A_{k,m}(1-A_{k,m})^{2}]\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}e_{k,i}^{2}(1)$$

$$+\frac{1}{n_{k}}E[A^{2}_{k,m}(1-A_{k,m})]\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}e_{k,i}^{2}(0)$$

$$-\frac{1}{n_{k}}E[A^{2}_{k,m}(1-A_{k,m})^{2}]\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}(e_{k,i}(1)-e_{k,i}(0))^{2}$$

$$=\frac{1}{n_{k}}E[A_{k,m}(1-A_{k,m})^{3}]\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}e_{k,i}^{2}(1)$$

$$+\frac{1}{n_{k}}E[A^{3}_{k,m}(1-A_{k,m})]\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}e_{k,i}^{2}(0)$$

$$+\frac{2}{n_{k}}E[A^{2}_{k,m}(1-A_{k,m})^{2}]\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}e_{k,i}(1)e_{k,i}(0)$$

$$=E\left[\frac{1}{n_{k}}\sum_{m=1}^{m_{k}}A^{3}_{k,m}(1-A_{k,m})^{3}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}\left(\frac{e_{k,i}(1)}{A_{k,m}}+\frac{e_{k,i}(0)}{1-A_{k,m}}\right)^{2}\right].$$

Therefore,

$$\liminf_{k\rightarrow\infty}E\left[\frac{1}{n_{k}}\sum_{m=1}^{m_{k}}A^{3}_{k,m}(1-A_{k,m})^{3}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}\left(\frac{e_{k,i}(1)}{A_{k,m}}+\frac{e_{k,i}(0)}{1-A_{k,m}}\right)^{2}\right]>0$$

is sufficient for $\liminf_{k\rightarrow\infty}f_{k}>0$ (even if condition (A.11) does not hold). Given our assumption that the supports of the cluster probabilities, $A_{k,m}$, are bounded away from zero and one (uniformly in $k$ and $m$), then

$$\liminf_{k\rightarrow\infty}E\left[\frac{1}{n_{k}}\sum_{m=1}^{m_{k}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}\left(\frac{e_{k,i}(1)}{A_{k,m}}+\frac{e_{k,i}(0)}{1-A_{k,m}}\right)^{2}\right]>0$$ (A.23)

is sufficient for $\liminf_{k\rightarrow\infty}f_{k}>0$. Assume that (A.23) holds, so $\liminf_{k\rightarrow\infty}f_{k}>0$. We now obtain,

$$E\Bigg{[}\Bigg{|}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i} W_{k,i}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})e_{k,i}(1)\Bigg{|}^{4}\,\Big{|}\,Q_{k,m}=1,A_{k,m}\Bigg{]}$$

$$=E\Bigg{[}(1-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})^{4}\Bigg{|}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}e_{k,i}(1)\Bigg{|}^{4}\,\Big{|}\,Q_{k,m}=1,A_{k,m}\Bigg{]}$$

$$\leq E\Bigg{[}\Bigg{|}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}e_{k,i}(1)\Bigg{|}^{4}\,\Big{|}\,Q_{k,m}=1,A_{k,m}\Bigg{]},$$

and

$$E\Bigg{[}\Bigg{|}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i} W_{k,i}e_{k,i}(1)\Bigg{|}^{4}\,\Big{|}\,Q_{k,m}=1,A_{k,m}\Bigg{]}$$

$$=E\Bigg{[}\Bigg{|}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}(R_{k,i}W_{k,i}-p_{k}A_{k,m})e_{k,i}(1)\Bigg{|}^{4}\,\Big{|}\,Q_{k,m}=1,A_{k,m}\Bigg{]}$$

$$=n_{k,m}E[(R_{k,i}W_{k,i}-p_{k}A_{k,m})^{4}|Q_{k,m}=1,A_{k,m}]$$

$$+3n_{k,m}(n_{k,m}-1)(E[(R_{k,i}W_{k,i}-p_{k}A_{k,m})^{2}|Q_{k,m}=1,A_{k,m}])^{2}.$$

The first equality holds because the terms $e_{k,i}(1)$ sum to zero within clusters. The second equality holds because, if $m_{k,i}=m_{k,j}=m$, with $i\neq j$, then $R_{k,i}W_{k,i}$ and $R_{k,i}W_{k,i}$ are independent conditional on $Q_{k,m}=1,A_{k,m}$, and $E[R_{k,i}W_{k,i}-p_{k}A_{k,m}|Q_{k,m}=1,A_{k,m}]=0$. Notice that

$$E[(R_{k,i}W_{k,i}-p_{k}A_{k,m})^{2}|Q_{k,m}=1,A_{k,m}]=p_{k}A_{k,m}(1-p_{k}A_{k,m})\leq p_{k},$$

which also implies $E[(R_{k,i}W_{k,i}-p_{k}A_{k,m})^{4}|Q_{k,m}=1,A_{k,m}]\leq p_{k}$. As a result,

$$\sum_{m=1}^{m_{k}}E\Bigg{[}\Bigg{|}\frac{1}{\sqrt{n_{k}p_{k}q_{k}}}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\} R_{k,i}W_{k,i}e_{k,i}(1)\Bigg{|}^{4}\Bigg{]}$$

$$\leq\frac{1}{n_{k}p_{k}q_{k}}+3\,\frac{\max_{m}n_{k,m}}{\min_{m}n_{k,m}}\frac{1}{m_{k}q_{k}}\rightarrow 0.$$

A.4 Derivations of the variance estimators

In this section, we derive the adjustments in the CCV variance. (We do this under the assumption that the $Z_{i}$ are independent. In our simulations we actually use a slightly different sampling scheme for the $Z_{i}$ where the average $\overline{Z}_{k,m}$ is identical and fixed in each cluster.) To derive the CCV variance of the least squares estimator, consider first a variance estimator of the form

$$\left(\sum_{i=1}^{n}V_{i}\right)^{2}.$$

We aim, however, to design an estimator based on a subsample consisting of units with $Z_{i}=1$, where $Z_{i}\in\{0,1\}$ is i.i.d. binary with $\Pr(Z_{i}=1)=p_{Z}$ and independent of $V_{i}$. First, notice that

$$E\left[\left(\sum_{i=1}^{n}V_{i}\right)^{2}\right]=\sum_{i=1}^{n}E[V_{i}^{2}]+2\sum_{i=1}^{n-1}\sum_{j=i+1}^{n}E[V_{i}V_{j}],$$

and

$$E\left[\left(\sum_{i=1}^{n}Z_{i}V_{i}\right)^{2}\right]=p_{Z}\sum_{i=1}^{n}E[V_{i}^{2}]+2p_{Z}^{2}\sum_{i=1}^{n-1}\sum_{j=i+1}^{n}E[V_{i}V_{j}].$$

Therefore,

$$E\left[\frac{1}{p_{Z}}\left(\sum_{i=1}^{n}Z_{i}V_{i}\right)^{2}\right]=\sum_{i=1}^{n}E[V_{i}^{2}]+2p_{Z}\sum_{i=1}^{n-1}\sum_{j=i+1}^{n}E[V_{i}V_{j}],$$

and

$$\frac{(1-p_{Z})}{p_{Z}^{2}}\left(E\left[\left(\sum_{i=1}^{n}Z_{i}V_{i}\right)^{2}\right]-p_{Z}\sum_{i=1}^{n}E[V_{i}^{2}]\right)=2(1-p_{Z})\sum_{i=1}^{n-1}\sum_{j=i+1}^{n}E[V_{i}V_{j}].$$

Adding the last two equations,

$$E\left[\left(\sum_{i=1}^{n}V_{i}\right)^{2}\right] =\frac{1}{p_{Z}^{2}}E\left[\left(\sum_{i=1}^{n}Z_{i}V_{i}\right)^{2}\right]-\frac{(1-p_{Z})}{p_{Z}}\sum_{i=1}^{n}E[V_{i}^{2}]$$

$$=\frac{1}{p_{Z}^{2}}E\left[\left(\sum_{i=1}^{n}Z_{i}V_{i}\right)^{2}\right]-\frac{(1-p_{Z})}{p_{Z}^{2}}\sum_{i=1}^{n}E[Z_{i}V_{i}^{2}].$$ (A.24)

The first term of the CCV variance estimator for least squares is based on the sample counterpart of the right-hand side of equation (A.24), with $1\{m_{k,i}=m\}R_{k,i}((W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k})\widehat{U}_{k,i}-(\widehat{\tau}_{k,m}-\widehat{\tau}_{k})\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k}(1-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k}))$ in the role of $V_{i}$.

To derive the CCV variance estimator for the fixed effect case, consider

$$\lambda_{k}=1-q_{k}\frac{(E[A_{k,m}(1-A_{k,m})])^{2}}{E[A_{k,m}^{2}(1-A_{k,m})^{2}]},$$

and let $f_{k}^{\rm CCV}=\lambda_{k}f_{k}^{\rm cluster}+(1-\lambda_{k})f_{k}^{\rm robust}$. This transformation is designed to reproduce the terms in $f_{k}$ with factor

$$\sum_{m=1}^{m_{k}}\frac{n^{2}_{k,m}}{n_{k}}(\tau_{k,m}-\tau_{k})^{2}.$$

These terms dominate $f_{k}$ as $k$ increases. It also reproduces several lower order terms.

Notice that

$$f_{k}^{\rm robust} =E[A_{k,m}(1-A_{k,m})^{2}]\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}e^{2}_{k,i}(1)+E[A^{2}_{k,m}(1-A_{k,m})]\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}e^{2}_{k,i}(0)$$

$$+\Big{(}E[A_{k,m}(1-A_{k,m})]-(5+p_{k})E[A^{2}_{k,m}(1-A_{k,m})^{2}]\Big{)}\sum_{m=1}^{m_{k}}\frac{n_{k,m}}{n_{k}}(\tau_{k,m}-\tau_{k})^{2}$$

$$+(2+p_{k})E[A^{2}_{k,m}(1-A_{k,m})^{2}]\sum_{m=1}^{m_{k}}\frac{n_{k,m}}{n_{k}}(\tau_{k,m}-\tau_{k})^{2}.$$

Then,

$$f_{k}^{\rm CCV}-f_{k} =(1-\lambda_{k})p_{k}E[A^{2}_{k,m}(1-A_{k,m})^{2}]\left(\sum_{m=1}^{m_{k}}\frac{n_{k,m}}{n_{k}}(\tau_{k,m}-\tau_{k})^{2}+\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}(e_{k,i}(1)-e_{k,i}(0))^{2}\right)$$

$$=p_{k}q_{k}(E[A_{k,m}(1-A_{k,m})])^{2}\left(\sum_{m=1}^{m_{k}}\frac{n_{k,m}}{n_{k}}(\tau_{k,m}-\tau_{k})^{2}+\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}(e_{k,i}(1)-e_{k,i}(0))^{2}\right).$$

For $\tilde{v}_{k}^{\rm CCV}=f_{k}^{\rm CCV}/(\mu_{k}(1-\mu_{k})-\sigma^{2}_{k})^{2}$, we obtain,

$$\tilde{v}_{k}^{\rm CCV}-\tilde{v}_{k}=p_{k}q_{k}\sum_{m=1}^{m_{k}}\frac{n_{k,m}}{n_{k}}(\tau_{k,m}-\tau_{k})^{2}+p_{k}q_{k}\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}(e_{k,i}(1)-e_{k,i}(0))^{2}.$$ (A.25)

The difference $\tilde{v}_{k}^{\rm CCV}-\tilde{v}_{k}$ is non-negative and of smaller order than $\tilde{v}_{k}$. Therefore, $\tilde{v}_{k}^{\rm CCV}/\tilde{v}_{k}\rightarrow 1$ (even if $\tilde{v}_{k}^{\rm CCV}-\tilde{v}_{k}$ is bounded away from zero). The first term on the right-hand side of (A.25) could be estimated to further correct the difference between the CCV estimator and the variance of $\widehat{\tau}_{k}^{\rm fixed}$.

A.5 Limit results

Let $X_{k,m}$ be an infinite array of random variables, with rows indexed by $k=1,2,\ldots$, and the columns of the $k$-th row indexed by $m=1,\ldots,m_{k}$. Let

$$S_{k}=\sum_{m=1}^{m_{k}}X_{k,m},$$

and $a_{k}=E[S_{k}]$.

A Weak Law of Large Numbers for Arrays: For each $k=1,2,\ldots$ , suppose that $X_{k,1},\ldots,$ $X_{k,m_{k}}$ are independent and have finite second moments. In addition, let $b_{k}$ be a sequence of positive constants such that

$$\frac{1}{b_{k}^{2}}\sum_{m=1}^{m_{k}}E[X_{k,m}^{2}]\longrightarrow 0.$$

Then,

$$\frac{S_{k}-a_{k}}{b_{k}}\stackrel{{\scriptstyle p}}{{\longrightarrow}}0.$$

Proof: By Chebyshev’s inequality, for any $\varepsilon>0$

$$\Pr\left(\left|\frac{S_{k}-a_{k}}{b_{k}}\right|>\varepsilon\right) \leq\frac{1}{b_{k}^{2}\varepsilon^{2}}\mbox{var}(S_{k})$$

$$=\frac{1}{b_{k}^{2}\varepsilon^{2}}\sum_{m=1}^{m_{k}}\mbox{var}(X_{k,m})$$

$$\leq\frac{1}{b_{k}^{2}\varepsilon^{2}}\sum_{m=1}^{m_{k}}E[X_{k,m}^{2}]\longrightarrow 0.$$

$\square$

A Central Limit Theorem for Arrays: For each $k=1,2,\ldots$ , suppose that $X_{k,1},\ldots,$ $X_{k,m_{k}}$ are independent, with zero means, $E[X_{k,m}]=0$, and finite variances, $\sigma^{2}_{k,m}=E[X^{2}_{k,m}]$, for $m=1,\ldots,m_{k}$. Let

$$s_{k}^{2}=\sum_{m=1}^{m_{k}}\sigma^{2}_{k,m}.$$

Assume also that Lyapounov’s condition holds,

$$\lim_{k\rightarrow\infty}\frac{1}{s_{k}^{2+\delta}}\sum_{m=1}^{m_{k}}E[|X_{k,m}|^{2+\delta}]=0,$$

for some $\delta>0$. Then,

$$\frac{S_{k}}{s_{k}}\stackrel{{\scriptstyle d}}{{\longrightarrow}}N(0,1).$$

Proof: billingsley, Chapter 27.

A.6 Intermediate calculations for Section A.2

The calculation of $v_{k}$ uses the following results.

$$E[(R_{k,i}W_{k,i}-p_{k}q_{k}\mu_{k})^{2}]=p_{k}q_{k}\mu_{k}(1-p_{k}q_{k}\mu_{k}),$$

$$E[(R_{k,i}(1-W_{k,i})-p_{k}q_{k}(1-\mu_{k}))^{2}]=p_{k}q_{k}(1-\mu_{k})(1-p_{k}q_{k}(1-\mu_{k})),$$

$$E[(R_{k,i}W_{k,i}-p_{k}q_{k}\mu_{k})(R_{k,i}(1-W_{k,i})-p_{k}q_{k}(1-\mu_{k}))]=-p_{k}^{2}q_{k}^{2}\mu_{k}(1-\mu_{k}),$$

$$E[R_{k,i}W_{k,i}R_{k,j}W_{k,j}|m_{k,i}=m_{k,j}]=E[p_{k}^{2}q_{k}A_{k,m}^{2}]=p_{k}^{2}q_{k}(\sigma_{k}^{2}+\mu_{k}^{2}),$$

and

$$E[(R_{k,i}W_{k,i}-p_{k}q_{k}\mu_{k})(R_{k,j}W_{k,j}-p_{k}q_{k}\mu_{k})|m_{k,i}=m_{k,j}] =p_{k}^{2}q_{k}(\sigma_{k}^{2}+\mu_{k}^{2})-(p_{k}q_{k}\mu_{k})^{2}$$

$$=p_{k}^{2}q_{k}(\sigma_{k}^{2}+(1-q_{k})\mu_{k}^{2}).$$

Similarly,

$$E[(R_{k,i}(1-W_{k,i}) -p_{k}q_{k}(1-\mu_{k}))(R_{k,j}(1-W_{k,j})-p_{k}q_{k}(1-\mu_{k}))|m_{k,i}=m_{k,j}]$$

$$=p_{k}^{2}q_{k}(\sigma_{k}^{2}+(1-q_{k})(1-\mu_{k})^{2}).$$

Notice also that

$$E[R_{k,i}W_{k,i}R_{k,j}(1-W_{k,j})|m_{k,i}=m_{k,j}] =E[p_{k}^{2}q_{k}A_{k,m}(1-A_{k,m})]$$

$$=p_{k}^{2}q_{k}(\mu_{k}(1-\mu_{k})-\sigma_{k}^{2}),$$

and

$$E[(R_{k,i}W_{k,i}-p_{k}q_{k}\mu_{k}) (R_{k,j}(1-W_{k,j})-p_{k}q_{k}(1-\mu_{k}))|m_{k,i}=m_{k,j}]$$

$$=p_{k}^{2}q_{k}(\mu_{k}(1-\mu_{k})-\sigma_{k}^{2})-p_{k}^{2}q_{k}^{2}\mu_{k}(1-\mu_{k})$$

$$=p_{k}^{2}q_{k}(\mu_{k}(1-\mu_{k})(1-q_{k})-\sigma_{k}^{2}).$$

  The following bounds are useful to prove Lyapunov’s condition.

$$E[|R_{k,i}W_{k,i}-p_{k}q_{k}\mu_{k}|^{3}] =(1-p_{k}q_{k}\mu_{k})^{3}p_{k}q_{k}\mu_{k}+(p_{k}q_{k}\mu_{k})^{3}(1-p_{k}q_{k}\mu_{k})$$

$$\leq c\,p_{k}q_{k}.$$

Let $Q_{k,m}$ be a binary indicator that takes value one if cluster $m$ of population $k$ is sampled.

$$E\big{[} |R_{k,i}W_{k,i}-p_{k}q_{k}\mu_{k}|^{2}|R_{k,j}W_{k,j}-p_{k}q_{k}\mu_{k}|\big{|}m_{k,i}=m_{k,j}=m\big{]}$$

$$=E\big{[}\big{(}(1-p_{k}q_{k}\mu_{k})^{2}p_{k}A_{k,m}+(p_{k}q_{k}\mu_{k})^{2}(1-p_{k}A_{k,m})\big{)}$$

$$\hskip 56.9055pt\times\big{(}(1-p_{k}q_{k}\mu_{k})p_{k}A_{k,m}+(p_{k}q_{k}\mu_{k})(1-p_{k}A_{k,m})\big{)}\big{|}m_{k,i}=m_{k,j}=m,Q_{k,m}=1\big{]}q_{k}$$

$$+E\big{[}\big{(}p_{k}q_{k}\mu_{k}\big{)}^{3}\big{|}m_{k,i}=m_{k,j}=m,Q_{k,m}=0\big{]}(1-q_{k})$$

$$\leq cp_{k}^{2}q_{k}.$$

$$E\big{[} |R_{k,i}W_{k,i}-p_{k}q_{k}\mu_{k}||R_{k,j}W_{k,j}-p_{k}q_{k}\mu_{k}||R_{k,t}W_{k,t}-p_{k}q_{k}\mu_{k}|\big{|}m_{k,i}=m_{k,j}=m_{k,t}=m\big{]}$$

$$=E\big{[}\big{(}(1-p_{k}q_{k}\mu_{k})p_{k}A_{k,m}+(p_{k}q_{k}\mu_{k})(1-p_{k}A_{k,m})\big{)}^{3}\big{|}m_{k,i}=m_{k,j}=m_{k,t}=m,Q_{k,m}=1\Big{]}q_{k}$$

$$+E\big{[}\big{(}p_{k}q_{k}\mu_{k}\big{)}^{3}\big{|}m_{k,i}=m_{k,j}=m_{k,t}=m,Q_{k,m}=1\Big{]}(1-q_{k})$$

$$\leq cp^{3}_{k}q_{k}.$$

 

Other useful intermediate calculations.

For the moments of treatment indicators, notice that $E[(W_{k,i}-\mu_{k})^{2}W_{k,i}]=\mu_{k}(1-\mu_{k})^{2}$, and $E[(W_{k,i}-\mu_{k})^{2}(1-W_{k,i})]=(1-\mu_{k})\mu_{k}^{2}$. In addition,

$$E[W_{k,i}W_{k,j}|m_{k,i}=m_{k,j}] =E[A_{k,m}^{2}]\quad\mbox{(for $m\in\{1,\ldots,m_{k}\}$)}$$

$$=\sigma_{k}^{2}+\mu_{k}^{2}.$$

Similarly, $E[(1-W_{k,i})(1-W_{k,j})|m_{k,i}=m_{k,j}]=\sigma_{k}^{2}+(1-\mu_{k})^{2}$. Therefore, $E[(W_{k,i}-\mu_{k})W_{k,j}|m_{k,i}=m_{k,j}]=\sigma^{2}_{k}$ and $E[(W_{k,i}-\mu_{k})(1-W_{k,j})|m_{k,i}=m_{k,j}]=-\sigma^{2}_{k}$. In addition,

$$E[(W_{k,i}-\mu_{k})(W_{k,j} -\mu_{k})W_{k,i}W_{k,j}|m_{k,i}=m_{k,j}]$$

$$=E[A^{2}_{k,m}](1-\mu_{k})^{2}\quad\mbox{(for $m\in\{1,\ldots,m_{k}\}$)}$$

$$=(\sigma_{k}^{2}+\mu_{k}^{2})(1-\mu_{k})^{2}.$$

Similarly,

$$E[(W_{k,i}-\mu_{k})(W_{k,j} -\mu_{k})(1-W_{k,i})(1-W_{k,j})|m_{k,i}=m_{k,j}]=(\sigma_{k}^{2}+(1-\mu_{k})^{2})\mu_{k}^{2},$$

and

$$E[(W_{k,i}-\mu_{k})(W_{k,j}-\mu_{k})W_{k,i}(1-W_{k,j})|m_{k,i}=m_{k,j}]=\mu_{k}(1-\mu_{k})(\sigma^{2}_{k}-\mu_{k}(1-\mu_{k})).$$

$\mbox{var}(R_{k,i}W_{k,i})=p_{k}q_{k}\mu_{k}(1-p_{k}q_{k}\mu_{k})$, $\mbox{var}(R_{k,i}(1-W_{k,i}))=p_{k}q_{k}(1-\mu_{k})(1-p_{k}q_{k}(1-\mu_{k}))$. Moreover,

$$\mbox{cov}(R_{k,i}W_{k,i},R_{k,i}(1-W_{k,i})) =E[R_{k,i}W_{k,i}R_{k,i}(1-W_{k,i})]-E[R_{k,i}W_{k,i}]E[R_{k,i}(1-W_{k,i})]$$

$$=-p_{k}^{2}q_{k}^{2}\mu_{k}(1-\mu_{k}).$$

Recall that $E[W_{k,i}W_{k,j}|m_{k,i}=m_{k,j}]=\sigma_{k}^{2}+\mu_{k}^{2}$. Therefore, $\mbox{cov}(W_{k,i},W_{k,j}|m_{k,i}=m_{k,j})=\sigma_{k}^{2}$. Also,

$$E[W_{k,i}(1-W_{k,j})|m_{k,i}=m_{k,j}]=\mu_{k}(1-\mu_{k})-\sigma_{k}^{2}.$$

$$E[R_{k,i}W_{k,i}R_{k,j}W_{k,j}|m_{k,i}=m_{k,j}] =E[R_{k,i}R_{k,j}|m_{k,i}=m_{k,j}]E[W_{k,i}W_{k,j}|m_{k,i}=m_{k,j}]$$

$$=p_{k}^{2}q_{k}(\sigma_{k}^{2}+\mu_{k}^{2}).$$

Similarly,

$$E[R_{k,i}(1-W_{k,i})R_{k,j}(1-W_{k,j})|m_{k,i}=m_{k,j}]=p_{k}^{2}q_{k}(\sigma_{k}^{2}+(1-\mu_{k})^{2}).$$

Therefore,

$$\mbox{cov}(R_{k,i}W_{k,i},R_{k,j}W_{k,j}|m_{k,i}=m_{k,j}) =p_{k}^{2}q_{k}(\sigma_{k}^{2}+\mu_{k}^{2})-p_{k}^{2}q_{k}^{2}\mu_{k}^{2}$$

$$=p_{k}^{2}q_{k}(\sigma_{k}^{2}+\mu_{k}^{2}(1-q_{k})),$$

and

$$\mbox{cov}(R_{k,i}(1-W_{k,i}),R_{k,j}(1-W_{k,j})|m_{k,i}=m_{k,j}) =p_{k}^{2}q_{k}(\sigma_{k}^{2}+(1-\mu_{k})^{2})-p^{2}_{k}q^{2}_{k}(1-\mu_{k})^{2}$$

$$=p_{k}^{2}q_{k}(\sigma_{k}^{2}+(1-\mu_{k})^{2}(1-q_{k})).$$

In addition,

$$\mbox{cov}(R_{k,i}W_{k,i},R_{k,j}(1-W_{k,j})|m_{k,i}=m_{k,j}) =E[R_{k,i}W_{k,i}R_{k,j}(1-W_{k,j})|m_{k,i}=m_{k,j}]$$

$$-E[R_{k,i}W_{k,i}|m_{k,i}=m_{k,j}]E[R_{k,j}(1-W_{k,j})|m_{k,i}=m_{k,j}]$$

$$=E[R_{k,i}R_{k,j}|m_{k,i}=m_{k,j}]E[W_{k,i}(1-W_{k,j})|m_{k,i}=m_{k,j}]$$

$$-E[R_{k,i}W_{k,i}|m_{k,i}=m_{k,j}]E[R_{k,j}(1-W_{k,j})|m_{k,i}=m_{k,j}]$$

$$=p_{k}^{2}q_{k}(\mu_{k}(1-\mu_{k})-\sigma_{k}^{2})-p_{k}^{2}q_{k}^{2}\mu_{k}(1-\mu_{k})$$

$$=p_{k}^{2}q_{k}(\mu_{k}(1-\mu_{k})(1-q_{k})-\sigma_{k}^{2}).$$

A.7 Intermediate calculations for Section A.3

 

$$E[R_{k,i}W_{k,i}(W_{k,i}-A_{k,m})|A_{k,m},Q_{k,m}=1,m_{k,i}=m]=p_{k}A_{k,m}(1-A_{k,m}).$$

This implies

$$E[R_{k,i}W_{k,i}(W_{k,i}-A_{k,m})|m_{k,i}=m]=p_{k}q_{k}E[A_{k,m}(1-A_{k,m})].$$

Therefore,

$$E\Bigg{[}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}(W_{k,i}-A_{k,m})\Bigg{]}=n_{k,m}p_{k}q_{k}E[A_{k,m}(1-A_{k,m})].$$

  For $n\geq 1$,

$$E\Bigg{[}\sum_{i=1}^{n_{k}} 1\{m_{k,i}=m\}R_{k,i}W_{k,i}(\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}-A_{k,m})\Big{|}\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$N$\kern-0.18004pt}}}_{k,m}=n\Bigg{]}$$

$$=\frac{1}{n}E\Bigg{[}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}-nA_{k,m}\Bigg{)}\Big{|}\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$N$\kern-0.18004pt}}}_{k,m}=n\Bigg{]}$$

$$=E[A_{k,m}(1-A_{k,m})].$$

Therefore,

$$E\Bigg{[}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}(\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}-A_{k,m})\Bigg{]} =E[A_{k,m}(1-A_{k,m})]\Pr(\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$N$\kern-0.18004pt}}}_{k,m}\geq 1)$$

$$=q_{k}E[A_{k,m}(1-A_{k,m})](1-(1-p_{k})^{n_{k,m}}).$$

 

For $n\geq 1$

$$E[R_{k,i}W_{k,i}(\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m} -A_{k,m})^{2}|m_{k,i}=m,\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$N$\kern-0.18004pt}}}_{k,m}=n,R_{k,i}=1]$$

$$\leq E[(\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}-A_{k,m})^{2}|m_{k,i}=m,\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$N$\kern-0.18004pt}}}_{k,m}=n,R_{k,i}=1]$$

$$\leq\frac{E[A_{k,m}(1-A_{k,m})]}{n}.$$

Because $\Pr(R_{k,i}=1|\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$N$\kern-0.18004pt}}}_{k,m}=n,m_{k,i}=m)=n/n_{k,m}$, we obtain

$$E[R_{k,i}W_{k,i}(\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}-A_{k,m})^{2}|m_{k,i}=m,\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$N$\kern-0.18004pt}}}_{k,m}=n]\leq\frac{E[A_{k,m}(1-A_{k,m})]}{n_{k,m}},$$

which implies

$$E[R_{k,i}W_{k,i}(\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}-A_{k,m})^{2}|m_{k,i}=m,\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$N$\kern-0.18004pt}}}_{k,m}\geq 1]\leq\frac{E[A_{k,m}(1-A_{k,m})]}{n_{k,m}}.$$

Therefore,

$$E[R_{k,i}W_{k,i} (\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}-A_{k,m})^{2}|m_{k,i}=m]$$

$$=E[R_{k,i}W_{k,i}(\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m}-A_{k,m})^{2}|m_{k,i}=m,\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$N$\kern-0.18004pt}}}_{k,m}\geq 1]\Pr(\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$N$\kern-0.18004pt}}}_{k,m}\geq 1|m_{k,i}=m)$$

$$\leq q_{k}\frac{E[A_{k,m}(1-A_{k,m})]}{n_{k,m}}.$$

 

Conditional on $\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$N$\kern-0.18004pt}}}_{k,m}=n$ and $A_{k,m}$, the variable $N_{k,m,1}$ has a binomial distribution with parameters $(n,A_{k,m})$. Then, using the formulas for the moments of a binomial distribution, we find that for any integer $n$, such that $1\leq n\leq n_{k,m}$,

$$E\Bigg{[}\Bigg{(}\sum_{i=1}^{n_{k}} 1\{m_{k,i}=m\}R_{k,i}W_{k,i}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})\Bigg{)}^{2}\Big{|}A_{k,m}=a,\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$N$\kern-0.18004pt}}}_{k,m}=n\Bigg{]}$$

$$=E[(N_{k,m,1}-N_{k,m,1}^{2}/n)^{2}|A_{k,m}=a,\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$N$\kern-0.18004pt}}}_{k,m}=n]$$

$$=n^{2}a^{2}(1-a)^{2}+na(1-a)(1-6a+6a^{2})+r_{1}(a)+r_{2}(a)/n,$$

where $|r_{1}(a)|$ and $|r_{2}(a)|$ are uniformly bounded in $a\in[0,1]$. Therefore,

$$E\Bigg{[}\Bigg{(}\sum_{i=1}^{n_{k}} 1\{m_{k,i}=m\}R_{k,i}W_{k,i}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})\Bigg{)}^{2}\Big{|}\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$N$\kern-0.18004pt}}}_{k,m}=n\Bigg{]}$$

$$=n^{2}E[A_{k,m}^{2}(1-A_{k,m})^{2}]+nE[A_{k,m}(1-A_{k,m})(1-6A_{k,m}+6A_{k,m}^{2})]$$

$$+E[r_{1}(A_{k,m})]+E[r_{2}(A_{k,m})]/n.$$

It follows that

$$E\Bigg{[}\sum_{m=1}^{m_{k}} (\tau_{k,m}-\tau_{k})^{2}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})\Bigg{)}^{2}\Bigg{]}$$

$$=\Bigg{(}\sum_{m=1}^{m_{k}}(\tau_{k,m}-\tau_{k})^{2}(n_{k,m}(n_{k,m}-1)p_{k}^{2}q_{k}+n_{k,m}p_{k}q_{k})\Bigg{)}E[A_{k,m}^{2}(1-A_{k,m})^{2}]$$

$$+\sum_{m=1}^{m_{k}}(\tau_{k,m}-\tau_{k})^{2}n_{k,m}p_{k}q_{k}E[A_{k,m}(1-A_{k,m})(1-6A_{k,m}(1-A_{k,m}))]+\mathcal{O}(m_{k}q_{k}).$$

Therefore,

$$\frac{1}{n_{k}p_{k}q_{k}}E\Bigg{[}\sum_{m=1}^{m_{k}} (\tau_{k,m}-\tau_{k})^{2}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i}(W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})\Bigg{)}^{2}\Bigg{]}$$

$$\longrightarrow (E[A_{k,m}(1-A_{k,m})]-(5+p_{k})E[A^{2}_{k,m}(1-A_{k,m})^{2}])\sum_{m=1}^{m_{k}}\frac{n_{k,m}}{n_{k}}(\tau_{k,m}-\tau_{k})^{2}$$

$$+p_{k}E[A^{2}_{k,m}(1-A_{k,m})^{2}]\sum_{m=1}^{m_{k}}\frac{n^{2}_{k,m}}{n_{k}}(\tau_{k,m}-\tau_{k})^{2}.$$

  Notice that,

$$E\Bigg{[}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i} (W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})\Bigg{)}^{4}\Big{|}A_{k,m}=a,\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$N$\kern-0.18004pt}}}_{k,m}=n\Bigg{]}$$

$$=E[(N_{k,m,1}(1-N_{k,m,1}/n))^{4}|A_{k,m}=a,\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$N$\kern-0.18004pt}}}_{k,m}=n]$$

$$\leq E[N^{4}_{k,m,1}|A_{k,m}=a,\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$N$\kern-0.18004pt}}}_{k,m}=n]$$

$$\leq n^{4},$$

Therefore,

$$E\Bigg{[}\Bigg{(}\sum_{i=1}^{n_{k}}1\{m_{k,i}=m\}R_{k,i}W_{k,i} (W_{k,i}-\hbox{\vbox{\hrule height=0.5pt\kern 1.1625pt\hbox{\kern-0.18004pt$W$\kern-0.18004pt}}}_{k,m})\Bigg{)}^{4}\Bigg{]}=n_{k,m}^{4}p_{k}^{4}q_{k}\left(1+\mathcal{O}\left(\frac{1}{p_{k}\min_{m}n_{k,m}}\right)\right),$$

uniformly in $m$.

 

Suppose $X_{k,m}=(Z_{k,m,1}+Z_{k,m,2})^{2}$. Let $X_{k,m,1}=Z_{k,m,1}^{2}$ and $X_{k,m,2}=Z_{k,m,2}^{2}$. Now suppose,

$$\sum_{m=1}^{m_{k}}E[X_{k,m,1}^{2}]\longrightarrow 0,$$

and

$$\sum_{m=1}^{m_{k}}E[X_{k,m,2}^{2}]\longrightarrow 0.$$

Using the binomial theorem and Hölder’s inequality, we obtain

$$\sum_{m=1}^{m_{k}}E[X_{k,m}^{2}] =\sum_{m=1}^{m_{k}}\sum_{p=0}^{4}c_{p}E[Z_{k,m,1}^{p}Z_{k,m,2}^{(4-p)}]$$

$$\leq c\sum_{m=1}^{m_{k}}\sum_{p=0}^{4}E[|Z_{k,m,1}|^{p}|Z_{k,m,2}|^{(4-p)}]$$

$$\leq c\sum_{m=1}^{m_{k}}\sum_{p=0}^{4}(E[X_{k,m,1}^{2}])^{p/4}(E[X_{k,m,2}^{2}])^{(4-p)/4}$$

$$\leq c\sum_{p=0}^{4}\left(\sum_{m=1}^{m_{k}}E[X_{k,m,1}^{2}]\right)^{p/4}\left(\sum_{m=1}^{m_{k}}E[X_{k,m,2}^{2}]\right)^{(4-p)/4}\longrightarrow 0.$$